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ABSTRACT 


Calculations  are  presented  that  are  elements  of  a  model  to  predict  the 
HR  radiation  emitted  by  high-altitude  rocket  plumes,  because  of  the  inter¬ 
action  of  the  plume  gases  with  the  ambient  atmosphere .  A  model  of  plume 
radiance  is  presented,  appropriate  to  high  altitudes  where  the  average  time 
between  plume  and  atmospheric  species  collisions  is  greater  than  the  radia¬ 
tive  lifetimes  of  the  relevant  excited  molecular  states.  The  three  parts  of 
the  calculation  -  flow  field  density,  molecular  collisional  excitation,  and 
molecular  nonequilibrium  radiation  -  are  identified  and  discussed. 

The  penetration  of  a  hypersonic  low-density  atmosphere  into  the  exhaust 
plume  of  a  steadily  firing  rocket  motor  is  studied  by  use  of  the  kinetic  theory 
of  gases.  The  flow  field  is  modeled  mathematically  by  employing  the  Hill- 
Draper  description  of  the  exhaust  jet  and  the  BGK  collision  integral  to  reduce 
the  problem  to  an  analytically  interpretable  and  computationally  tractable 
form.  The  distribution  function  of  the  ambient  gas  is  bimodal,  with  the 
unscattered  fraction  of  the  molecules  peaked  about  the  uniform  freestream 
velocity,  and  the  scattered  gas  carried  along  with  the  macroscopic  exhaust 
gas  velocity. 

Classical  rotational  excitation  probabilities  for  the  collisions  of  HF, 

HCl,  and  OH  with  three  atomic  species  have  been  calculated  in  three  dimen¬ 
sions  by  use  of  a  Monte  Carlo  procedure.  Results  are  presented  for  the 
rotational  excitation  of  both  rigid  diatomic  rotors  and  those  coupled  via  a 
classical  harmonic  oscillator.  Hard-sphere  intermolecular  potentials  ai-e 
employed  to  obtain  results  applicable  to  high-energy  ( ~1  eV)  collisions. 

The  results  of  these  two  methods  of  calculation  indicate  that  for  the  center- 
of-mass  collisional  energy  range  of  1  to  10  eV,  which  dominates  the  problem 
of  interest,  activation  of  the  vibrational  modes  is  quite  inefficient. 

Once  excited  to  high  rotational  states,  these  diatomic  molecules  can 
radiate  at  infrared  wavelengths .  The  power  radiated  into  a  frequency  band 
via  pure  rotational  radiation  by  a  diatomic  molecule  with  a  permanent  electric 
dipole  moment  has  been  calculated  both  quantum  mechanically  and  classically. 
The  results  define  the  region  of  validity  for  the  classical  approximation.  It 
is  concluded  that  for  the  range  of  parameters  of  interest  to  the  plume  inter¬ 
ference  problem  the  classical  approximation  is  adequate.  In  order  to  be  con¬ 
crete,  the  results  of  the  theory  are  applied  to  the  particular  case  of  rotational 
emission  from  the  hydrogen  flouride  molecule. 
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FAR  FIELD  PLUME  RADIANCE  MODEL 


1..  INTRODUCTION 

Plume  interference  can  be  roughly  defined  as  the  degradation  of  the  signal-to- 
noise  ratio  of  a  sensor  due  to  the  radiation  emitted  by  the  exhaust  plume  of  a  rocket 
engine  either  aboard  the  same  vehicle  as  the  sensor  or  aboard  some  other,  nearby 
vehicle.  If  the  engine  is  aboard  the  same  vehicle  as  the  sensor,  the  situation  is  one 
of  self-interference;  if  it  is  aboard  a  nearby  vehicle,  it  is  one  of  mutual  interference. 
The  emphasis  of  the  work  reported  here  is  on  self-interference,  although  much,  if 
not  all,  of  it  is  applicable  to  mutual  interference  as  well. 

In  the  summer  of  1971,  a  panel  with  Dr.  Hans  Wolfhard  of  IDA  as  chairman 
performed  a  preliminary  study  of  the  plume  interference  problem!1^  In  particular, 
this  panel  divided  plume  interference  into  tv-o  basic  parts:  the  near-field  radiation 
and  the  far-field  radiation.  The  near-field  was  defined  as  the  radiation  originating 
from  IR-active  molecular  species  emerging  from  the  rocket  nozzle  in  an  excited 
state.  The  far-field  was  defined  as  the  radiation  from  plume  species  excited  into 
IR-active  states  because  of  collisions  with  atmospheric  species.  Thus,  the  division 
of  the  problem  is  based  on  the  molecular  excitation  mechanism . 

On  the  basis  of  their  initial  examination  of  the  plume  interference  problem, 
the  Wolfhard  panel  recommended  that  it  be  studied  in  more  detail.  The  first  efforts 

to  study  the  problem  undertaken  by  ABMDA  were  directed  towards  the  specific  case 

(2  3) 

of  the  proposed  HIT  interceptor  system  .  ’ 

This  report  deals  with  the  far-field  plume  interference  problem,  and  is  the 
Final  Technical  Report  of  the  far-field  modeling  efforts  performed  at  Aerodyne 
Research,  Inc,  during  the  ABMDA  Plume  Interference  Program.  The  work  done  is 
a  natural  extension  of  far-field  plume  radiance  calculations  made  earlier  under  the 

(4) 

HIT  Plume  Interference  Program. 
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The  calculations  performed  to  predict  the  far-field  plume  radiance  of  the  HIT 
vehicle  were  for  a  small,  short-pulsed  solid  propellant  engine  fired  in  the  upper 
atmosphere,  typically  above  200  km.  The  molecular  radiating  mode  identified  as 
the  dominant  source  of  radiation  in  the  far  field  was  pure  rotation  of  water.  The 
details  of  the  various  calculations  are  discussed  in  the  HIT  Program  Final  Technical 
Report^  as  well  as  in  separate  publications^.5,6,7^ 

These  HIT  program  calculations  were  performed  over  a  relatively  short  period 
of  time  in  order  to  calculate  the  signal-to-noise  ratio  degradation  as  a  function  of 
time  after  engine  shutdown,  and  the  results  of  the  calculations  were  used  as  system 
design  limits  and  criteria.  The  work  reported  below,  and  performed  under  the  Plume 
Interference  Program,  was  done  under  different  circumstances,  however.  The  aim 
of  the  far-field  modeling  under  this  later  program  was  to  broaden  the  predictive 
capability  of  the  initial  modeling  efforts  by  considering  different  aerodynamic  regimes, 
and  by  including  IR-active  molecular  modes  other  than  the  water  rotation  considered 
for  HIT.  To  achieve  this  end  a  number  of  calculations  were  undertaken.  (In  addition, 
data  analyses  and  plume  radiance  predictions  were  carried  into  execution,  and  are 
reported  under  separate  cover  in  Volumes  II  and  III). 

Section  3  of  this  report  discusses  a  calculation  of  the  penetration  of  a  hypersonic 
low-density  atmosphere  into  the  exhaust  plume  of  a  steadily  firing  rocket  motor  using 
the  kinetic  theory  of  gases.  Section  4  concerns  the  far-field  molecular  excitation 
process;  classical  rotational  excitation  probabilities  for  collisions  of  HF,  HC1,  and 
OH  with  three  atomic  species  have  been  calculated  in  three  dimensions.  Section  5 
concerns  the  nonequilibrium  radiation  from  the  excited  molecular  states;  the  power 
radiated  into  a  frequency  band  via  pure  rotational  radiation  by  a  diatomic  molecule 
with  a  permanent  electric  dipole  moment  has  been  calculated  both  quantum  mechanically 
and  classically.  The  calculations  presented  in  these  three  sections,  along  with  the 
overall  plume  radiance  formulation  discussed  in  Section  2,  are  the  far-field  plume 
radiance  model  elements  performed  during  this  program. 

Before  proceeding  to  the  details,  it  is  perhaps  worthwhile  to  outline  the  present 
status  and  extent  of  the  far-field  models.  Thus  far,  there  are  basically  two  aerodynamic 
interaction  configurations  which  can  be  modeled.  The  first  configuration  is  the  decay 
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of  short-pulsed  engine  plumes  operating  at  high  altitudes.  The  engines  aboard  the 
proposed  HIT  vehicle  are  of  such  a  type.  The  second  is  the  interaction  region  between 
the  relatively  dense  core  of  a  continuously  firing  engine  and  the  high-altitude  atmo¬ 
sphere.  In  each  case,  the  engine  size  appropriate  to  the  model  has  up  to  several 
thousand  pounds  thrust  at  altitudes  of  400  to  500  km . 

In  each  of  these  two  aerodynamic  regimes,  a  number  of  IR-active  plume  species 
have  been  included  in  the  overall  radiation  model.  The  high  collision-energy  excita¬ 
tion  of,  and  the  subsequent  non  equilibrium  radiation  from  the  pure  rotational  motion 
of  H20  and  the  diatomic  hydride  molecules  HF,  HC1  and  OH  have  so  far  been  con¬ 
sidered.  Concentration  has  been  centered  on  these  pure  rotational  emitters,  since 
rotational  radiation  forms  a  quasicontinuous  spectral  background  across  the  wave¬ 
length  regions  of  interns'  in  the  plume  interference  problem,  and  thus,  cannot  be 
eliminated  by  judicious  *  ivelength  filtering. 

While  considering  uie  scope  of  the  present  far-field  plume  radiance  models,  it 
is  worth  noting  that  all  t  ie  calculations  mentioned  have  begun  from  basic  principles. 
There  is  verv  little,  if  any,  data  presently  available  shedding  light  on  the  far-field 
problem.  Thus,  the  modeler  bears  a  heavy  burden,  and  it  behooves  him  to  base 
his  models  as  firmly  as  possible  on  accepted  physical  principles. 


r>  v  A  *  * 
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2.  PLUME  RADIANCE  FORMULATION 


The  overall  objective  of  the  far-field  modeling  efforts  is  to  be  able  to  calculate 
plume  radiance  for  a  variety  of  interceptor  configurations.  In  order  to  handle 
this  objective,  the  overall  radiance  problem  has  been  broken  down  into  several  parts, 
as  was  mentioned  in  the  introduction.  This  section  will  more  clearly  define  this 
division  into  parts,  and,  in  particular,  will  describe  the  overall  radiance  formulation 
necessary  to  combine  these  parts. 

Generally  speaking,  the  sequence  of  events  leading  to  far-field  radiation 
from  a  rocket  exhaust  plume  is  as  follows.  The  rocket  engine  burns  for  some  time, 

At,  and  an  expanding  mass  of  gas  is  emitted  from  the  nozzle.  As  this  gas  flows  away 
from  the  nozzle,  the  plume  molecular  species  begin  to  collide  with  atmospheric 
species.  During  these  collisions,  plume  species  are  excited  into  IR  active  states, 
and  they  consequently  begin  to  radiate.  As  they  continue  to  travel  away  from  the 
rocket  nozzle,  the  plume  species  continue  to  undergo  collisions  that  can  excite  them 
(or  possibly  deexcite  them  if  they  have  not  had  time  to  radiate  away  the  excitation 
energy  imparted  in  a  previous  collision).  In  addition,  these  collisions  bring  the 
plume  molecules  into  equilibrium  with  the  ambient  atmosphere,  so  that,  asymptotically, 
in  time,  the  plume  becomes  a  cloud  of  gas  at  rest  with  respect  to  the  atmosphere 
and  diffuses  away  Therefore,  the  three  individual  parts  of  the  far-field  radiance 
problem  are  the.  calculations  of  the  plume  and  atmosphere  gas  densities,  the  excita¬ 
tion  rate  of  relevant  modes  of  plume  species,  and  the  power  radiated  by  these  excited 
modes.  The  p  ume  radiance  model,  then  combines  these  separate  calculations  to 
calculate  plume  radiance. 

The  overall  radiance  model  to  be  used  here  is  similar  to  the  one  evolved  earlier 
for  the  specific  case  of  the  HIT  intercept  vehicled  The  following  discussion  parallels 
the  discussion  of  that  modeHo  be  found  in  Ref.  (8). 

The  plume  is  an  extended  source  of  radiation.  The  power  on  the  detector  from 
such  an  extended  source  is  given  simply  by 
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w'h_re  £1  f  is  the  solid  angular  field  of  view  of  the  detector,  A„  is  the  area  of  the 
I  K 

receiver  optics,  and  N  is  the  plume  radiance  (emitted  into  the  operating  band  of  the 

The  source  radiance  can  be  written 

dS.  I  (S.)  n  (S.  ;•  r,  t)  (2-2) 

where  I  (S^  is  the  power  radiated  into  the  band  by  a  water  molecule  in  the  STh  state, 
n  (S^  ;  r  ,  t)  is  the  density  of  such  molecules  at  point  r  in  space  and  at  time  t,  the 
line  integral  dr  is  along  the  sensor  line  of  sight,  and  the  integral  dS.  is  over  the 
manifold  of  accessible  molecular  states  S. .  Although  correct,  this  equation  is  not 

i 

very  useful  as  written.  For  the  conditions  under  which  many  interceptor  systems 
operate,  considerable  simplification  is  possible. 

First,  at  high  altitudes  (e.g. ,  above  200  km)  the  atmosphere  is  sufficiently 
tenuous  that  the  average  time  between  collisions  of  plume  and  atmospheric  species 
becomes  longer  than  the  radiative  lifetimes  of  the  relevant  LWIR  active  excited  states, 
and  most  excited  molecules  completely  deexcite  via  radiative  emission  between 
collisions.  This  limit  of  short  radiative  lifetimes  compared  to  mean  collision  times 
is  the  exact  opposite  limit  of  the  more  usually  encountered  one  of  thermal  equilibrium. 

Second,  for  a  narrow  field-of-view  detector,  the  plume  and  atmospheric  gases  may 
be  considered  uniform  across  the  field  of  view.  This  means  that  the  power  radiated 
by  an  ensemble  of  molecules  in  some  elemental  volume  is  given  by 

dP  =  EN  (E*?  Wj,  aig)  n  (E*  ;  r,  t  )  (2-3) 

where  n  is  the  creation  rate  of  molecules  in  an  excited  state  E*,  and  EN  is  the  total 

$ 

energy  ultimately  radiated  by  the  state  E  into  the  frequency  band  of  interest,  uj^  to 


2 

detector;  units:  watts/cm  -sr). 

OO 

N  -  ~tr\  dr  / 
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Thus,  to  find  the  plume  radiance  under  these  conditions,  one  need  only  add  up  the 
contributions  of  all  the  excited  states  and  integrate  along  the  observational  line  of 
sight.  For  the  case  of  rotational  excitation  and  radiation,  of  primary  interest  here 
as  explained  in  the  introduction,  the  plume  radiance,  N,  can  be  written  as 

OO  OO 

N  =  “  f  dr  f  dER  EN  (ER  ;  cur  cu2)  n  (ER ;  7,  t)  (2-4) 

0  (/ 

where  FTT  (ER;  cuj,  cj2)  is  the  total  uj^  —  band  energy  emitted  by  a  molecule 

in  the  initial  rotational  energy’  state  E„,  n  is  the  rate  of  creation  of  such  molecules 

,  K 

at  point  r  and  at  time  t  per  unit  energy,  the  integral  dr  is  still  a  line  integral  along 
the  line  of  sight,  and  the  dER  integral  is  over  all  initial  rotational  energies  of  the 
molecule. 

n  (E„  ;  r  ,  t),  the  rate  of  creation  of  molecules  in  the  state  E,,  per  unit  volume 

K  K 

per  unit  energy  (E„),  is  related  to  the  differential  excitation  cross  section,  the 

collision  relative  velocity,  and  the  exhaust  and  ambient  gas  densities.  That  is,  n 

is  the  number  of  collisions  between  exhaust  molecules  and  ambient  atmospheric 

particles  per  unit  energy  per  unit  time  producing  the  appropriate  final  rotational 

state,  and  is  given  by 


^  (Er  >  f, 


'/A  a/' 


(V  )  f 
a  e 


where  vg  is  the  ambient  particle  velocity,  v  is  the  plume  molecule  velocity, 
da  /d  E„  is  the  differential  rotational  excitation  cross  section,  f,  (v  )  is  the 
ambient  atmospheric -phase  space  distribution  function  (units:'  cm-3  (cm/sec)  ), 
and  fg  (v  )  is  the  exhaust -phase  space  d'stribution  function  for  the  appropriate 
plume  constituent.  This  expression  is  the  usual  collisional  rate  constant  expression 
for  an  ambipolar  process.  For  the  moment,  the  two  distribution  functions  will 
remain  unspecified  to  maintain  complete  generality,  but  it  is  worthwhile  to  note 
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that  in  virtually  every  conceivable  case  of  interest,  the  phase  space  distribution 
function  is  separable  into  a  spatial  distribution  function  (i.e. ,  number  density)  and 
a  velocity  distribution  function.  For  example,  for  a  uniform  ambient  atmosphexe 
in  thermal  equilibrium  at  temperature  ,  fg  (v  )  would  simply  become 


f 
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oo 
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2ir  kT 


exp 


m  v 
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2  kT 


(2-6) 


where  n  is  a  constant  number  density  (in  the  context  of  Eq.  (2-6)),  and  the  rest  of 
the  expression  is  the  usual  Maxwellian  velocity  distribution. 

Combining  the  expression  for  h,  Eq.  (2-5),  with  the  expression  for  the  plume 
radiance  N,  Eq.  (2-4),  yields  the  following: 
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If  the  two  distributions  are  separable  into  two  distinct  pieces. 
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then 


N  = 


/ 

LOS 


dr  -4 ST  J  dER  J  d\  j  d3ve 
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va  "  ve  6e  <ve)  ga  (va) 


(2-9) 


-*■  _3 

where  g  (v)  is  a  velocity  distribution  function  (units:  (velocity)  ),  ng  is  the 
number  denstiy  of  the  relevant  plume  species,  n  is  the  ambient  number  density, 
and  ng  and  na  (expressed  as  ne  a  (r,  t))  are  both  functions  of  time  and  position. 

This  expression  for  the  plume  radiance,  N,  along  with  the  detailed  calculations 
of  EN  (E„;  cu. ,  co0),  dcr/9  E  ,  n  and  n  ,  form  the  basis  of  the  present  far-field 
plume  background  interference  analysis.  The  results  of  the  separate  calculations  of 
n^,  ng,  acr/d  E^,  and  EN  can  be  used  in  an  expression  like  Eq.  (2-7)  to  calculate 
a  plume  radiance  for  a  sensor  line  of  sight.  This  formulation  01  the  plume  radiance 
problem  s  appropriate  for  calculating  the  power  falling  on  some  given  detector.  If 
more  detailed  information  is  desired,  plume  radiance  contours  and  a  somewhat 
different  formulation  are  required.  In  the  present  context,  however,  it  should  be 
emphasized  that  the  individual  calculations  to  be  described  in  the  following  sections 
of  this  report  are  independent  of  the  details  of  the  exact  radiance  formulation,  and 
can  be  used  in  more  detailed  plume  radiance  calculations,  or  in  radiance  formula¬ 
tions  that  relax  the  short  radiative  lifetime  or  spatially  homogeneous  assumptions 
discussed  above. 
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3.  ATMOSPHERIC  PENETRATION  OF  EXHAUST  PLUMES  UNDER 
RARIFIED  FLOW  CONDITIONS 

3. 1  Introduction 

The  flow  pattern  set  up  by  the  interaction  of  an  exhaust  plume  generated  by  the 

continuous  firing  of  a  rocket  motor  and  a  tenuous  atmosphere  has  been  intens;vely 

studied  in  recent  years.  Most  of  this  work  has  focused  on  those  domains  of  motor 

thrust  and  operating  altitude  permitting  this  interaction  to  be  described  on  the  basis 

of  continuum  fluid  mechanics.  At  altitudes  sufficiently  high  for  the  ambient  mean 

free  path  to  become  larger  than  the  continuum  interaction  length  scale,  the  analysis 

of  such  a  flow'  must  be  based  on  the  kinetic  theory  of  gases.  Under  these  conditions, 

the  length  scale  associated  w'ith  the  atmospheric  penetration  of  the  exhaust  plume  is 

\  A/Kn,  where  A  is  the  nozzle  exit  area,  and  Kn  is  the  Knudsen  number  based  on 

V^and  the  mean  free  path,  based  on  average  exit  conditions  (Muntz,  Hamel,  and 
(9) 

Maguire  1970);  '  This  scale  is,  typically,  much  larger  than  the  dimensions  of  the 
vehicle  carrying  the  motor,  since  Kn  «  1.  The  presence  of  the  vehicle  itself  may 
thus  be  ignored  in  analyzing  phenomena  on  this  scale. 

When  the  penetration  scale  is  small  compared  with  the  ambient  mean  free  path, 
the  response  of  the  atmosphere  to  the  presence  of  the  plume  is  decoupled  from  the 
process  by  which  the  atmosphere  modifies  the  vacuum  plume  expansion.  This 
situation  is  appropriate  to  moderate-sized  engines,  such  as  those  employed  for 
attitude  control  or  maneuvering  purposes  in  the  upper  atmosphere.  The  spatial 
separation  of  these  two  phenomena  has  been  exploited  by  Brook  and  Hamel  (1972)^*^ 
in  their  study  of  a  spherical  source  interacting  with  a  stationary  background  gas. 
Baum  (1973)^^  considered  the  expansion  of  a  transient  exhaust  plume  into  a 
rarefied  atmosphere  by  invoking  equivalent  assumptions  about  the  scales  associated 
with  the  interactions  occurring  in  that  problem. 


PRECEDING  PAGE  BLANK-NOT 
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The  present  work  is  concerned  with  the  penetration  of  a  hypersonic  ambient 

strain  into  a  steady-state  exhaust  plume  under  the  conditions  stated  above.  The 

exhaust  gas  properties  are  computed  from  a  simple  model  proposed  by  Hill  and 
(12) 

Draper  (1966),  '  of  the  continuum  expansion  of  a  jet  into  a  vacuum..  This  expanding 

jet  serves  as  a  diffuse  scatterer  of  ambient  molecules.  The  scattering  is  represented 
mathematically  by  the  Krook  collision  integral.  This  formulation  leads  to  an 
analytically  interpretable  and  computationally  tractable  integral  equation  for  the 
ambient  density  distribution..  The  molecular  distribution  function  is  described 
analytically  in  terms  of  the  density.  The  distribution  function  is  bimodal,  with  the 
unscattered  molecules  peaked  about  the  uniform  free-stream  velocity,  and  the  scattered 
gas  convected  with  the  local  macroscopic  exhaust  velocity.  An  asymptotic,  analytical 
solution  to  the  integral  equation  is  obtained.  This  result  is  used,  in  conjunction  with 
direct  numerical  solutions,  to  obtain  illustrative  flow  patterns  for  both  axially  symme¬ 
tric  and  nonsymmetric  configurations. 


3.2  Mathematical  Formulation 

The  formulation  of  the  model  will  now  be  considered  in  detail.  The  penetration 
scale  is  assumed  to  be  sufficiently  small  compared  with  the  ambient  mean  free  path 
for  collisions  between  ambient  molecules  to  be  ignored..  The  kinetic  equation  for 
the  ambient,  molecular  distribution  function  f(x,  v)  is  a  modified  version  of  that  given 
by  Gross  and  Krook  (1956). 
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Here  n  (x)  is  the  ambient  number  density,  and  ne>  u^,  Tg  are,  respectively,  the 
exhaust  gas  number  density,  macroscopic  velocity,  and  temperature.  The  maximum 
velocity  of  the  exhaust  gas  is  u^,  while  the  cross  section  cr,  entering  the  expression 
for  the  collision  frequency  ce,  is  an  effective  momentum  transfer  cross  section. 

The  relation  between  the  momentum  transfer  between  species  computed  from  Eq.  (3-1), 
and  that  given  by  Gross  and  Krook  (1P56)^^  is  explained  in  Baum  (1973)!1  ^  The 
mean  relative  speed  is  taken  to  be  that  between  u^  and  the  ambient  wind  ,  because 
the  exhaust  and  ambient  speed  ratios ,  Wg  and  ,  are  large  under  the  conditions  of 
interest. 

The  exhaust  properties  are  assumed  in  a  modified  version  of  the  form  devised 
by  Hill  and  Draper  (1966)(12^  The  adopted  forms  are  equivalent  to  the  Hill-Draper 
model  far  from  the  engine,  but  take  on  the  prescribed  nozzle  exit  values  at  x  =  0. 

The  relevant  equations  are: 


BA  ..  exp 
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The  parameters  B  and  A 


are  given  in  terms  of  exit  properties  by 
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The  nozzle  exit  number  density  n^,  speed  u^,  temperature  Tq,  Mach  number 
and  the  ratio  of  specific,  heats  y  are  determined  by  the  engine  parameters.  The 
exhaust  properties  are  expressed  in  a  spherical,  polar  coordinate  system,  centered 
about  the  unit  thrust  vector  t.  The  polar  angle  is  9 ,  and  the  azimuthal  angle  is 
measured  from  the  plane  containing  the  wind  and  thrust  vectors.  The  wind  vector 
makes  an  angle  '/'with  the  thrust  vector.  Thus  (see  Fig.  3-1) 

A  A  A  A 

r  (t  cos 6  +  j  sin0cos<£  +  k  sinBsixup )  =  r  r 

A  a 

u  (t  cos\p  +  j  smy) 

A  A  A  .  r  .  A 

Vv  =  V  (cos£t  +  sin^cos^j  +  sincsin^k) 

All  quantities  with  a  caret  superscript  are  unit  vectors.  The  system  of  Eqs.  (3-1) 
through  (3-3)  must  be  supplemented  by  the  boundary  condition  that,  far  from  the 
plume,  the  flow  is  uniform. 


L  f(x,  v) 

r  —*■  oc 


Here  n  ,  n  ,  and  T  are,  respectively,  the  ambient  number  density,  velocity, 

oo  oo  oo 

and  temperature  far  from  the  plume. 

The  representation  of  the  exhaust  gas  displayed  in  Eq.  (3-2)  and  (3-3)  describes 
the  expansion  into  vacuum  of  a  continuum  jet  plume.  The  motion  is  isentropic  and 
is  energetic.  The  density  distribution  (and,  hence,  each  displayed  fluid  quantity) 
satisfies  the  inviscid  equations  of  motion  asymptotically,  far  from  the  nozzle  exit. 
The  angular  dependence  of  the  asymptotic  density  distribution  is  chosen  as  an 
empirical  curve  fit  to  numerical  methods  of  characteristic  solutions  for  a  variety  of 
engines.  The  expressions  for  the  parameters  B  and  A  ensure  global  conservation 
of  mass,  momentum,  and  energy. 
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The  continuum  description  of  the  exhaust  gas  ultimately  breaks  down  as  the  plume 

(14) 

expands.  The  rarefaction  process  has  been  studied  by  Grundy  (1969);  '  who  shows 

that,  except  at  large  angles  Q ,  with  respect  to  the  thrust  axis,  the  density  and  velocity 
are  given  by  the  continuum  solution.  The  radial  component  of  the  pressure  tensor,  on 
the  other  hand,  does  not  continue  to  decrease  isentropically,  but  levels  off  at  a  finite, 
"frozen"  asymptote.  The  Mach  number  corresponding  to  the  freezing  temperature 
is  typically  very  large,  on  the  order  of  10-20  for  most  angles  Q  .  Under  these 
circumstances  ,  the  interaction  between  ambient  and  exhaust  gases  is  nearly  indepen¬ 
dent  of  the  local  plume  temperature,  as  shown  by  the  assymptotic  solution  attained 
in  Section  3.4. 

The  next  step  in  the  calculation  is  the  construction  of  the  fundamental  integral 
equation  for  n(x).  Following  Anderson  (1965);  ’  Eq.  (3-1)  may  be  written  in 

characteristic  form  as 


This  describe,0,  the  evolution  of  f,  at  any  point  r,  a  distance  s  along  a  ray  through 
the  field  point  x,  m  the  direction  v  (Fig.  3-2).  Now  integrate  Eq.  (3-  5)  from  minus 
infinity  to  the  field  point  s=o,  applying  the  boundary  condition  represented  by  Eq.  (3-4). 
The  result  is 
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Note  that  the  V  dependence 

$  depends  only  on  the  angle  variables  C  and  ^7  . 
Eq.  (3-6),  by  using  the  definition  of  n  in  terms 


The  argument  of  ng,  V  and  Te  in  Eq.  (3-6)  is  x  +  vs’ 
is  now  displayed  explicitly ,  since 
The  integral  equation  follows  from 
of  the  distribution  function  f . 
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Thus 
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cos/3  =  cosijcos^  +  sin^sinsj/cost) 


Again,  all  exhaust  properties  in  Eq.  (3-7)  are  functions  of  x+vs'.  The  functions 
Hn  (p  ,q)  are  defined  by: 
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Anderson  and  Macomber  (1964)'  '  contains  a  detailed  analysis  of  this  function, 

together  with  tables  for  n  =  1,  2,  3.  Eq.  (3-7)  determines  the  spatial  density 
distribution  of  the  atmosphere  as  it  penetrates  the  plume  core,  while  Eq.  (3-6) 
expresses  the  distribution  function  in  terms  of  n.  In  this  form,  the  equations  are 
too  complicated  to  be  useful.  However,  by  utilizing  the  fact  that,  in  many  cases  of 
interest,  both  Wg  and  are  large,  Eqs.  (3-6)  and  (3-7)  can  be  reduced  to  a 
tractable  form. 
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3.3  The  Hypersonic  Limit 


Consider  the  first  term  in  the  integral  Eq.  (3-7).  For  large  W  ,  the  expo¬ 
nential  factor  suggests  that  the  integrand,  considered  as  a  function  of  is  strongly 
2 

peaked  about  cos/3  =  1.  Examination  of  the  Hn  functions  reveals  that  Hn,  con¬ 
sidered  as  a  function  of  p,  is  exponentially  small  for  large  negative  p  and  propor- 
n — 2 

tional  to  p  for  large  positive  p.  Since  the  dependence  of  H  on  W  for  large 

^  n— 2 

positive  p  is  weak,  part  from  a  purely  multiplicative  factor  W  ,  the  dominant 
contribution  to  the  £  integration  should  come  from  the  neighborhood  of  cos/3  =  +1. 

To  proceed  formally  from  these  ideas  it  is  more  convenient  to  do  the  £  integral 
first,  and  employ  the  variable  x=cos£  .  The  integral  under  consideration  assumes 
the  form : 
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Applying  the  method  of  steepest  descent  to  (3-8),  the  saddle  point  is  at: 
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Both  saddle  points  are  real,  but  the  contribution  from  the  one  associated  with  the 
minus  sign  is  exponentially  smaller  than  that  associated  with  the  positive  value  of 
coscr,  and  is  henceforth  ignored.  The  integral  1  is  now  x’eadily  evaluated  as: 
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Using  the  result  (3-9),  the  integration  over  rj  may  be  readily  performed.  The 
required  integral  is  then: 


o 


Examination  of  this  expression  when  W  sin<j  »  1  (i.e. ,  the  undisturbed  flow  is 
neither  parallel  nor  anti-parallel  to  the  thrust  vector)  shows  that  the  dominant 
contribution  comes  from  the  saddle  point  at  =  o.  The  apparent  saddle  point  at 
t)  =  ir  is  spurious,  with  the  result  that: 


T  =  t  d  -  0  rj  -  0)  (3-10) 
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A  separate  calculation  when  sin  $  is  0(1)  (but  Woc»l)  yields  the  same  result  for  J. 

Physically,  Eq.  (3-10)  states  that  the  overwhelming  majority  of  the  unscattered 
ambient  molecules  are  aligned  with  the  unifoim  wind  upstream  of  the  plume.  A 
similar  analysis  (the  details  of  which  are  presented  in  Appendix  A)  may  be  applied  to 
the  homogeneous  term  in  Eq.  (3-7).  The  hypersonic  limit  of  Eq.  (3-7)  then  become o; 
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The  inhomogeneous  term  in  Eq.  (3-11)  contains  the  variable  =  \  —  t  (x+as), 

where  a  is  a  unit  vector  in  the  direction  of  the  undisturbed  wind.  This  quantity  is 
directly  proportional  to  the  number  of  collisions  experienced  by  a  molecule  traveling 
in  a  straight  line  from  infimcy  to  the  point  in  question  in  the  direction  of  the  undisturbed 
wind.  Written  explicitly,  it  has  the  form: 
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As  y  -*«  ,  the  number  of  collisions  vanishes  and  the  inhomogeneous  term  in 
Eq.  (3-11)  approaches  unity.  As  y  —  o,  the  number  of  collisions  experienced  becomes 
very  large  and  the  inhomogeneous  term  vanishes  exponentially.  The  number  of 
collisions  also  becomes  large  when  <p  =  o  and  6  =y  ,  since  the  molecule  has  then  had 
to  traverse  the  near  singularity  in  exhaust  gas  density  at  the  origin.  Physically,  the 
ray  on  the  opposite  side  of  the  nozzle  exit  from  the  wind  is  shielded  from  the  ambient 
molecules  for  finite  distances  r  from  the  exit.  The  nature  of  q^  near  the  singular  ray 
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corresponding  to  cos^g*  =  1  may  be  analyzed  as  follows:  The  behavior  is  dominated 
by  the  singularity  at  cos /S  *  =  1,  as  y  approaches  infinity  or  Kn  vanishes.;  Thus,  the 
quantities  cos/z  and  cosS  in  Eq.  (3-13)  may  be  replaced  by  their  values  at  cos =  1, 
i.  e. 

cos/z  =  -  cos\p,  1  +  t  <  0 
cos/z  =  cos 1  +  t  >0 
cosS  =  -1  ,  1  +  t  <0 
cosS  =  1  ,  1  +  t  >  0 

Upon  inserting  these  values  into  Eq.  (3-13)  the  integral  for  q^  may  be  evaluated 
explicitly  as: 
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This  expression  is  exact  along  the  singular  ray  cos/3*  =  1,  where  it  reduces  to: 


Along  the  singular  ray,  the  hypersonic  approximation  breaks  down  for  large  y,  since 
the  boundary  condition  at  infinity  requires  to  vanish  there.  Actually,  unscattered 
molecules  arriving  from  directions  other  than  that  of  the  freestream  flow  reach  the 
ray  9=^,  <p  =  0  for  large  y,  permitting  the  boundary  condition  to  be  satisfied.  Along 
any  adjacent  ray  corresponding  to  a  fixed  value  of  cos/3*  such  that  cos/3*  ^  1,  when 

n  2 

(Kn/y)  is  sufficiently  small  compared  withl  -  cos  /3*,  Eq.  (3-14)  reduces  to: 


Thus,  the  boundary  condition  at  infinity  is  recovered,  but  with  increasing  slowness  as 

the  singular  ray  is  approached.  Note  that  the  quantity  y  \l  -  cos2/3*  is  proportional 

to  the  distance  in  the  plane  perpendicular  to  the  ray  /3*  -  1  from  the  ray  to  the  point 

in  question.  The  ’'shadow”  region  is  then  a  cylindrical  tube  surrounding  the  singular 

_  2 

ray.  Since  the  inhomogeneous  time  (W^  ,  q^  )  is  roughly  an  exponential 

in  (“  /W^ ),  the  width  of  the  shadow  region  is  given  (in  units  Vba  Kn  *)  by  the 

2  -1 

1-cos  /3*)  in  Eq.  (3-16). 
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Now  consider  the  homogenous  term  in  Eq.-  (3-li).  The  variable  qg  =(m/kl'e)r* 
has  the  explicit  form: 
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The  quantity  qg  is  proportional  to  the  number  of  collisions  experienced  by  a  molecule 

initially  penetrating  to  a  distance  y'  along  a  given  radial  ray  as  it  travels  outward 

along  the  ray  to  the  field  point  at  y.  Note  that  q  /W  (y1)  becomes  very  large  (0  (Kn  *) 
o  e  e 

fore  <  60  and  realistic  values  of  and  u^  /u  )  at  y'  =  0  for  fixed  y.  Since  the 

function  H3  (we.Qe)  is  roughly  exponential  in  (-  q  A'e),  the  probability  that  those  few 

molecules  deposited  upstream  of  th  <it  can  penetrate  to  the  other  side  is  very  small. 


Moreover,  the  contribution  to 


ity  at  any  point  from  molecules  scattered  radially 


inwardtowards  the  nozzle  is  exj.  .lly  smaller  in  W  than  that  from  molecules 
scattered  radially  outward  away  i  .he  nozzle.:  For  these  reasons,  the  lower  limit 
of  integration  in  the  homogenous  term  of  Eq.  (3-11)  is  the  origin.  Similarly,  the 
overwhelming  preponderance  of  outward  scattering  events  means  that  there  is  no 
contribution  to  the  number  density  at  y  from  points  y*  >  y  on  that  ray.  Further 
details  are  given  in  Appendix  A.,  Finally,  the  functions  Hn  may  be  simplified  b\ 
applying  tl e  method  of  steepest  descent  once,  again.  The  result,  valid  for  p  >>  1  and 

Z  =  q/)?  ^  finite,  is  the  first  term  in  Anderson  and  Macomber’s  (1964)^^  asymptotic- 
expansion: 
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The  quantity  p  Z  in  the  exponential  is  equal  to  q^  /W  in  the  inhomogeneous  term  of 

Eq.  (3-11),  and  to  q  /W  in  the  homogeneous  term.  These  quantities  are  independent 

of  VV  and  W  ,  respectively.  Thus,  the  quantity  W  and  W  enter  the  H  functions 
oo  e  J  M  °°2  e-2  n 

only  through  the  quantity  Z,  which  is  formally  )  or  0(\Ve  )  for  fixed  v. 

These  terms  are  retained  in  Eq.  (3-18)  to  allow  y  to  vary  from  a  moderate  multiple 

of  Kn  to  values>>  1.  The  principal  effect  of  this  variation  is  to  shift  the  location  of 

the  saddle  point  tQ  for  radius  of  y  <  Although  the  shift  is  not  large,  the  presence 

of  the  quantity  t  in  the  exponent  of  Eq.  (3-18)  can  change  the  computed  value  of 

n/n^  at  a  given  location  by  as  much  as  a  factor  of  two. 

The  physical  picture  represented  by  Eq.  (3-11)  may  be  summarized  as  follows: 
Molecules  entering  in  a  beam  parallel  to  the  undisturbed  flow  are  scattered  out  of  the 
beam  by  the  exhaust  gas,  which  is  unaffected  by  the  process.  The  scattered  molecules 
are  thus  deposited  along  each  radial  ray  in  accordance  with  the  inhomogeneous  term 
in  Eq.  (3-11).  The  molecules  are  then  caught  up  in  the  radial  expansion  and  swept 
away  from  the  nozzle,  undergoing  additional  collisions  as  they  proceed  radially  out¬ 
ward.  The  speed  ratio  of  both  the  incident  beam  and  the  scattered  gas  is  so  high  that 
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the  thermal  spread  in  molecular  velocities  may  be  ignored,  in  comparison  with  the 
velocities  considered.  The  exhaust  gas  density  is  so  high  in  the  vicinity  of  the  nozzle 
exit  that  no  molecules  can  penetrate  this  region.  As  the  radial  distance  from  the 
nozzle  exit  increases,  the  penetration  process  becomes  more  efficient.  Ultimately, 
the  exhaust  gas  density  is  sufficiently  low  for  the  ambient  molecules  to  pass  through 
the  plume  undisturbed,  and  the  ambient  density  is  recovered.  The  phenomena  outlined 
above  are  illustrated  schematically  in  Fig.  3-3. 

The  remainder  of  this  paper  deals  with  the  solution  of  Eq.  (3-11),  and  the 
computation  of  two  examples.  The  following  section  contains  an  asymptotic,  analytic 
solution  of  this  equation,  valid  far  from  the  nozzle  exit.  Section  3.5  discusses  the 
computational  procedures  employed  in  the  direct  numerical  solution  of  Eq.  (3-11)  and 
in  the  evaluation  of  the  asymptotic  result  of  Section  3.6.  The  analytical  and  numerical 
methods  explained  in  these  two  sections  are  applied  to  an  axially  symmetric  flow 
W  =  0)  and  a  nonsymmetric/'side  blowing"  flow  (\|/  =  90°)  in  Section  3. 6. 


3.4  The  Asymptotic  Solution 

When  the  radial  coordinate  y  is  0(1),  the  solution  to  Eq.  (3-11)  may  be  readily 
obtained.;  Under  these  circumstances,  ue(y’)  approaches  u^,  while  q  /W  becomes: 
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Thus,  q  /W  is  0(11,  which  means  that  the  quantity  Z  =  q  /W  2  in  the  function  H„ 

C  C  0  0  O 

may  be  set  equai  to  zero,  reducing  Hg  to  an  exponential.  Eq,  (3-11)  then  reduces  to: 
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Schematic  of  Nonlocal  Interaction 
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Equation  (3-20)  is  readily  solved  by  differentiating  with  respect  to  y  to  obtain  a  first- 
order  differential  equation  for  the  quantity  jn/n^  -  1/W^  H^(y)J  J(y).  Examination 
of  Eq.  (3-20).  as  y  approaches  zero,  shows  that  this  quantity  should  vanish  at  the 
origin.  Thus: 
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The  function  H4(y)  can  be  simplified  in  a  similar  manner  away  from  the  singular  ray. 
The  function  q^  then  has  the  form: 
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Then,  provided  that  the  numerical  value  of  .0,  on  the  ray  under  consideration  is  not 

3  2 

too  large,  q^  /W^  may  be  set  equal  to  zero  in  the  function (l/W^H^  (W^  q^  ), 

thereby  reducing  it  to  an  exponential  in  (-  il/y).  The  quadrature  on  the  righthand  side 

of  Eq.  (3-21)  may  then  be  evaluated  in  terms  of  the  exponential  integral  function 

E^fZ).  yielding  the  result: 
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The  first  term  in  the  solution  (3-23)  could  also  form  ally  be  replaced  by  exp  j-,Q./yj  . 
to  the  present  order  of  accuracy.  However,  the  analysis  in  this  section  can  be 
interpreted  as  an  attempt  to  find  an  approximate  expression  for  the  scattered  gas 
contribution  to  the  number  density,  as  represented  b\  the  homogeneous  term  in 
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Eq.  (3-11).  The  second  term  in  the  solution  (3-23)  is  the  desired  result,  as  is  the 
righthand  side  of  Eq.  (3-21).  However,  the  other  expression  is  too  complicated  to 
study  analytically,  and  sufficiently  expensive  in  computer  time  (relative  to  the  direct 
numerical  solution  of  Eq.  (3-11))  to  warrant  its  use  without  the  further  approximation 
represented  by  Eq.  (3-23). 

Now  consider  the  behavior  of  the  solution  as  a  function  of  y  and  6 .  For  small 
y  (assuming  for  the  moment  that  the  asymptotic  solution  is  at  least  qualitatively 
correct)  the  function  E^{Sl/y)  has  the  asymptotic  expansion 


E, 


-  m 
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-  ii/y 


i  +o 


(f) 


The  scattered  gas  contribution  becomes  very  small  as  the  nozzle  exit  is  approached 
because  few  ambient  molecules  are  able  to  penetrate  the  dense  region  of  the  plume. 
For  large  angles  9,  the  factor  exp  j-  Aj*(l-cos0)2j  (representing  the  angular 
dependence  of  the  exhaust  gas  density  distribution)  cuts  off  the  scattered  gas  contri¬ 
bution  because  there  are  so  few  exhaust  gas  molecules  to  serve  as  scatterers.  For 
large  y,  E ^  (£l/y)  may  be  approximated  by: 


Thus,  the  scattered  gas  density  rises  to  a  peak,  then  decays  much  more  slowly  with 
increasing  y.  The  decay  is  caused  by  the  increasing  volume  in  the  plume  available 
to  molecules  scattered  much  nearer  to  the  nozzle  ana  then  convected  radially  outward. 
The  additional  local  scattering  far  from  the  nozzle  contributes  weakly,  accounting 
for  the  logarithmic  factor  in  the  scattered  result.  However,  the  exhaust  gas  density 
ultimately  decreases  to  the  point  where  the  unscattered  contribution  dominates, 
permitting  the  uniform  boundary  condition  at  infinity  to  be  recovered. 
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3. 5  The  Numerical  Computation  Procedure 

The  numerical  computation  of  n/nM  requires  the  evaluation  of  the  inhomogeneous 
term  in  Eq.  (3-11)  or  Eq.  (3-23),  together  with  a  procedure  for  calculating  the 
homogeneous  term  in  the  integral  equation.  The  evaluation  of  the  asymptotic  approxi¬ 
mation  to  the  Hn  functions  is  straightforward,  as  is  the  computation  of  the  function  E2, 
as  a  function  of  its  argument.  Hence  attention  is  confined  to  the  computation  of  q^  , 
which  determines  the  value  of  the  inhomogeneous  term;  and  to  the  evaluation  of  the 
homogeneous  term  in  Eq.  (3-11). 

The  computation  of  the  integral  defining  q^  will  be  considered  first.  Note  that 
the  integrand  is  not  singular  in  the  entire  range  of  integration,  since  its  denominator  is 
equal  to  or  greater  than  one  throughout  this  range.  However,  as  is  shown  in  Section  3.3, 
the  limit  of  q^  as  y  and  cos/3*  1  is  non-uniform;  that  is,  the  value  of  the  limit 
depends  on  the  order  m  which  these  two  limiting  processes  take  place.  Numerically, 
this  non-uniform ity  means  that  one  must  be  very  careful  in  evaluating  this  integral 
when  both  y  is  0(1)  or  larger  and  cos  0*  is  very  close  to  1.  For  such  values  of  y  and 
cos/3  *,  most  of  the  contributions  to  the  integral  come  from  points  very  close  to  t  =  -1, 
and  the  accuracy  with  which  cosS  and  cos^-  can  be  evaluated  at  these  points  is  limited 
by  the  arithmetical  precision  of  a  computer. 

A  close  examination  of  the  variation  of  the  integrand  J  of  Eq.  3-13  with  respect  to  t 

shows  that  the  integrand  varies  smoothly  in  the  entire  range  of  integration,  so  long  as 
*  * 

cos/3  is  not  close  to  1.  For  cos/3  near  1,  the  integrand  has  very  sharp  variations 

* 

near  t  =  -1,  and  a  peak  value  in  that  neighborhood.  For  cos 0  =  1,  the  peak  of  the 

integrand  is  exactly  at  t  =  -1,  and  is  discontinuous. 

In  viev'  of  the  above  described  peculiarities  at  the  integrand  J,  the  iollowing 
method  was  devised  to  compute  the  integral.  The  integral,  I,  is  divided  into  two 
parts,  IL  and  IR: 
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1  =  *L  +  *R 
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where  ^  and  are  two  small  positive  numbers,  is  chosen  such  that  the  omitted 
part  of  the  integral,  from  -l-€j  to  -1  divided  by  1^,  is  equal  to  the  specified  fractional 
tolerance  (TOL)  for  the  computation  of  the  inhomogeneous  and  homogeneous  terms. 
is  chosen  in  a  similar  fashion.  It  can  he  shown  that,  for  cos/3*  =  1,  and  are 
given  by: 
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Equations  (3-24)  are,  in  fact,  used  to  compute  and  e2  for  any  value  of  cos/3*..  It 
can  be  shown  that  these  values  of  and  are  smaller  than  the  required  values,  thus 
resulting  in  smaller  values  for  the  omitted  parts  of  the  integral.. 
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1^  is  computed  in  steps  At  of  the  integration  variable  t,  starting  from  t  =  -l-€ 
and  advancing  to  the  left  towards  t  =-°°  .  The  magnitude  of  each  step  At  is  determined, 
such  that  the  integrand  changes  by  about  one  order  of  magnitude  in  the  range  of  the 
step.  The  contribution  to  1^  from  each  step  is  completed  through  a  six-point  Gaussian 
integration.  1^  is  computed  in  a  similar,  steplike  fashion,  starting  from  t  =  -1 
and  advancing  to  the  right  towards  the  origin.  1^  and  1^  are  computed  concurrently, 
so  that  the  total  integral  I  is  progressively  built  up  from  both  sides  of  the  point  t  =  -1. 
The  computation  of  1^  is  terminated  when  the  most  recent  contribution  to  1^,  divided 
by  the  accumulated  I,  is  less  than  the  fractional  tolerance,  TOL,  of  the  computation. 
The  computation  of  1^  is  terminated  in  a  similar  fashion,  or  when  the  origin  is  reached. 

Now  consider  the  calculation  of  the  homogeneous  term  in  the  integral  equation. 

For  computational  purposes,  it  is  more  convenient  to  work  with  the  radial  coordinate 
r*  ,  defined  by: 


r*  =  — I —  =  (Kn)  *y 
V  BA 


(3-25) 


The  initial  growth  of  the  homogeneous  term  in  this  variable  is  slow  and  well-behaved, 
taking  place  over  the  range  10  <  r*  <  100  along  most  radial  rays.  Letting  the  dimen¬ 
sionless,  inhomogeneous  and  homogeneous  quantities  in  Eq.  (3-11)  be  n.*  and  n^*, 
respectively,  this  equation  may  be  written  in  the  form: 


n  *(r')  =  n.*  (r ’)  +  n^*  (r') 


(3-27) 
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In  the  "exact"  definition  of  the  integral  equation,  r2*  is  equal  to  zero.  However, 
in  the  numerical  evaluation  of  the  integral  of  this  equation  r^*  can  be  much  larger 
than  zero  as  will  be  shown  later. 

Now,  consider  the  following  mesh  of  points: 


y*  *  r  ♦  p  ♦ 

rl  •  x2  *  3  * 


r»+  r*+  Y*+ 

•  m-2’  m-l*  m’  m+1 


r*  ,  =  r*  +  h 
m+1  m  m 


(3-28) 


Note  that,  in  general,  a  variable  mesh  size  is  assumed.  Assume  that  the  values  of 
n*  are  known  at  all  points  up  to  and  including  the  point  r^  .  The  value  of  n*  at  tne 
next  point,  r^+1,  is  computed  by  the  following  procedure. 

The  integral  of  equation  (3-26)  is  divided  uTo  two  parts  so  that: 
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(3-31) 


Ij  is  evaluated  by  the  parabolic  integration  formula  (3-B10)  of  Appendix  B. 
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I 


] 


Thus: 


I-.  =  - 


■xr - 7t - t - :  K  (r*  ,,  r*  1)n*(r*  .) 

6h  ,  (h  t-  h  ,)  '  m-1  m+1'  v  m-1 

m-1  '  m  m-1 


h  (3h  +  h  ) 

m  m-1  m 


6  h 


m-1 


K  (r*  ,  r*  )  n*  (r*  ) 
'  m  m+1  '  m 


h  (3h 
m  '  m- 


,  +  2  h  ) 
1  m 


6(h 


m 


+  h  t) 

m-i 


K  (r*  •  -i » 

'  m+1’ 


rm+l> 


[ni  (rm+l> 


n*  (r  .) 
h  '  m+17 


(3-32) 


Substituting  Ij  from  Eq.  (3-32)  in  Eq.  (3-29),  and  solving  the  resulting  equation  for 
n^(r*m+i)i  one  obtains: 


II  +  I, 


n£<* 


h'  m+1 


)  - 


(3h 

m'  m- 


,  +  2  h  )/6(h  +  h 
1  nr  '  m  m 


-1*  K  *rm+l’  rrn+l) 


i 

l 

I 

I 


where  1^  is  et  ual  to  I ^ ,  as  given  by  Eq.  (3-32),  with  n£(r^1+i)  set  to  zero. 

As  can  be  .'.een  from  Eqs.  (3-31)  and  (3-32),  the  righthand  side  of  Eq  (3-33) 
involves: 

(i)  The  values  of  n*  for  all  points  up  to  and  including  the  point  r^,  which  are 
known;  and 

(ii)  The  value  of  n’t:  at  the  point  r^+j  which  can  i  •  evaluated  separately. 

Thus,  the  righthand  side  of  Eq.  (3-33)  involves  quantities  that  are  known,  or  can  be 
evaluated  independently  of  n£  (r^+1)’’  therefore,  Eq.  (3-33,  is  an  explicit  expression 
for  the  unknown  nj^  (rJ^+1)- 


I 
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Ig  is  computed  by  successive  applications  of  the  parabolic  integration  formula 
(3-B11)  of  Appendix  B  to  three  points  at  a  time,  starting  with  the  point  r^  and 
advancing  to  the  left  towards  the  point  r*.  The  reason  for  computing  I2  backwards, 
in  this  way,  is  that  the  kernel  K(r',  r^+1)  decreases  sharply  as  r'  becomes  smaller 
and  smaller  than  r*  ,  and  thus  the  contributions  to  I„  decrease  sharply.  This  means 
that  the  computation  of  I2  can  be  terminated  when  a  desired  accuracy  is  achieved, 
much  before  reaching  the  initial  point  r*.  The  computation  of  I2  is  terminated  when 
the  most  recent  contribution  to  I2,  divided  by  the  sum  of  I2  (accumulated)  and  I'  ,  is 
smaller  than  a  specified  fractional  tolerance  (TOL),  which  is  also  the  tolerance 
within  which  n|  is  computed.  If  m  is  even,  and  the  computation  of  I2  must  be  extended 
all  the  way  to  r*,  then  the  last  contribution  to  I2  is  from  two  points,  rather  than 
three,  namely,  points  r*  and  r*  ,  and  formula  (3-B9)  of  Appendix  B  is  applied  instead 
of  formula  (3-B11). 

The  procedure  described  above  for  the  computation  of  n£  requires  that  n£  be 
known  at  the  first  three  points,  rj ,  r£  and  r^.  The  value  of  n£  at  r|  is  the  initial 
condition  of  the  problem,  and  when  r|  is  taken  equal  to  zero,  one  has  the  ’’exact 
initial  condition"  of  the  problem,  which  is  nj^(O)  =  0..  A  close  examination  of  the 
situation  reveals  that,  as  one  travels  along  a  radial  ray  from  the  origin  outward,  n* 
decays  rapidly  from  1  to  a  very  small  value,  then  increases  back  to  1,  while  nj* 
increases  slowly  from  0  to  some  maximum  value,  then  decays  back  to  zero.  Thus, 
in  the  numerical  evolution  of  n£,  a  great  deal  of  computer  time  is  saved,  without  any 
loss  in  accuracy,  by  starting  the  solution  of  some  r*  larger  than  zero,  with  the 
approximate  initial  condition  n£  =  0  at  r* .  r*  is  selected  so  that  it  is  in  the  region 
where  n’t  has  its  minimum  plateau.  The  rapidly  decaying  kernel  annihilates  the 
errors  associated  with  this  approximate  initial  condition.  Thus: 


n*(r|)  =  0 


(3-34) 
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The  value  of  n£  at  r*  is  computed  by  evaluating  lj  in  fcq.  (3-30)  through  the 
trapezoidal  rule  (Ig  in  this  case  is,  of  course,  equal  to  zero).  As  for  the  value  of 
nh  r3’  **  *s  comPutec*  by  evaluating  both  1^  and  Ig  in  Eqs.  (3-30)  and  (3-31)  by 
the  trapezoidal  rule.  Thus: 


n*  <r|) 


h2/2  [K(r|,  r|)  n*(rx)  +  K(r*,  r|)  n^(r*)] 
1  -(h1/2)K(r|,  r|) 


(3-35) 


nh(r3)  = 


1/2  [hx  K(r*,  r*)  n*(r*)  +  <h 


1+V 


K(r2’ 


r*)  n*(r* )  + 


h2K(r3’  r3>  ni  <r3> 


1  -(h„/2)K(r*  r*) 


(3-36) 
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If  the  integrating  step,  hm,  from  r^  to  r^+  ^  is  too  large,  then  the  denominator 
in  Eq.  (3-33)  can  be  small,  or  even  negative.  This  must  not  be  allowed  since 
instabilities  will  result.  Now,  it  can  be  shown  that: 


/ 


\ 


#  *  *  \  ~  /  *  * 

K  r  .  i ,  r  xi  I  <  K  (  r  ,  r  ] 

'  m+ 1  m+ 1 j  \m  m ) 


(3-37) 


and,  therefore,  if  the  integrating  step  is  chosen  such  that:. 


h  = 
m 


a 


K  (  r*  ,  r*  ^ 
'  m  m  j 


(3-38) 


then  the  denominator  of  Eq.  ^3-33)  is  roughly  equal  to  one  to  0.5a.  Thus,  for  suffi¬ 
ciently  small  values  of  the  parameter  a,  instabilities  can  be  avoided..  The  correct 
value  of  g 1  for  a  given  radial  ray  cannot  be  defined  a  priori,  but  one  anticipates  that  a 
value  between  0. 1  and  1  should  do.  Indeed,  in  the  applications  discussed  in  the  next 
section,  the  value  a  =  0.5  was  employed.  This  ambiguity  in  the  right  value  of  a 
notwithstanding,  Eqs.  (3-37)  and  (3-38)  show  that  hm  increases  as  rj^  increases,  a 
very  desirable  result.  Another  criterion  for  the  selection  of  hm  is  derived  by  fitting 
a  parabola  to  the  values  of  n*  at  the  three  points,  r*  r*  ,,  and  r*  ,  then  extrapo- 
lating  to  a  point  r^+  where  the  fractional  change  in  nj^  is  a  small  number  /?,  on 
the  order  of  0.01  to  0. 1.  The  value  of  /?  actually  used  is  0.05.  The  resulting  value 
of  h^  is  compared  to  that  given  by  Eq.  (3-38),  and  the  smaller  of  the  two  is  selected. 
For  the  two  initial  integrations  with  the  trapezoidal  rule,  a  sufficiently  small  integrat¬ 
ing  step  is  taken  readily  satisfying  the  stability  considerations  discussed  above. 

The  computer  program  is  in  FORTRAN  IV,  and  has  been  run  successfully  on  the 
UNIVAC  1108  system.  Before  final  production  runs  were  made,  the  accuracy,  stabil¬ 
ity,  and  efficiency  of  the  program  were  investigated  through  several  test  runs,  by  use 
of  a  trial  integral  equation  whose  exact  solution  was  known. 


•  t'J  /v*i  V&ih*fT?l3  »*H-&,  >  fl-?*? 
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The  procedure  for  calculating  n.  works  very  well  for  any  value  cf  r  ,  so  long  as 
cos  [i  is  not  equal  to  1.  For  cos  /3  =1,  the  computer  program  is  only  accurate 
for  moderate  values  of  r.  However,  the  solution  is  not  physically  meaningful  for 
r*  sufficiently  large  for  y  to  be  0(1)  on  this  ray.  The  procedure  for  calculating  n* 
is  stable  and  accurate.  It  is  efficient  for  values  of  0  equal  to  or  less  than  about  45° 
if  an  accuracy  of  0. 1%  (TOL  =  0.001)  is  required.  By  accepting  an  error  of  0.5% 

(TOL  =  0.005),  the  program  can  be  utilized  out  to  angles  6  of  60°.  For  larger  values 
of  6,  the  value  of  the  parameter  a,  required  for  an  accurate  solution,  is  so  small  that 
the  integrating  steps  become  very  small,  and  a  great  deal  of  computer  time  is  required 
to  obtain  the  solution.  For  such  values  of  9,  n^  can  be  computed  much  more  efficiently 
from  the  asymptotic  expression  (3-23),  valid  quite  close  to  the  origin  for  such  large 
angles.  A  computer  program  with  the  asymptotic  expressions  for  n*  has  also  been 
developed,  and  works  quite  well. 

For  an  accuracy  of  one  part  in  a  thousand  (TOL  =  0.001),  it  takes,  on  the  average, 

about  30  sec  of  1108  central  processing  unit  (CPU)  time  to  compute  n*  and  n*  along  a 

radial  ray,  all  the  way  from  the  origin  to  the  value  of  r*  where  the  ambient  condition 

is  recovered.  The  CPU  time  required  for  the  calculation  of  n*  is  about  three  times 

that  for  nf. 
i 

3.6  Numerical  Results 

The  theoretical  and  numerical  analyses  described  in  the  previous  sections  have 
been  applied  to  the  calculation  of  an  axially  symmetric  flow  ("aligned  flow")  $  =  0: 
and  to  a  "non-aligned"  case,  tp  =  90°.  The  values  of  the  remaining  parameters  needed 
to  uniquely  determine  the  solution  for  n/nw  were  the  same  for  both  flows,  and  are 
listed  in  Table  3-1. 

TABLE  3-1.  Parameters  Employed  in  Computed  Results 


Kn 

\v 

oo 

UM/U~ 

y 

M 

o 

m/m 

Kf3 

6 

0.3 

1.3 

4 

0.8 
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The  integral  equation  was  solved  directly,  to  an  accuracy  of  0. 1%,  on  rays  spaced  at 
five-degree  integrals  for  0  <45°,  and  with  a  reduced  accuracy  of  0.5%  for  values  of  0 
between  45°  and  60°.  Equation  (3-23)  was  ased  to  obtain  results  for  angles  0  greater 
than  60°. 

First  consider  the  aligned  flow:  Figure  (3-4)  illustrates  the  relative  magnitude 
of  the  scattere  i  and  unscattered  contributions  to  the  density  on  the  0  =  30°  ray.  Near 
the  nozzle  exit.,  both  contributions  are  small  because  most  of  the  molecules  in  the 
incident,  ambient  stream  are  deposited  before  reaching  the  30°  ray.  As  y  increases, 
more  molecules  are  deposited  on  the  ray,  and  then  convected  outward  by  the  exhaust. 
Since  the  scattered  gas  at  a  given  y  station  is  an  accumulation  of  molecules  scattered 
earlier,  this  juildup  is  quite  rapid,  and  gives  rise  to  the  overshoot  in  density  depicted 
in  the  figure.  The  combination  of  the  geometric  relieving  effect  and  reduced  scattering 
discussed  in  Section  3.4  then  reduces  the  scattered  gas  contribution  very  slowly  (note 
the  logarithmic  scale)  to  zero.  The  unscattered  gas  density  rises  monotonic  ally  with 
increasing  y  to  recover  the  boundary  condition  at  infinity. 

The  abcve  picture  holds  qualitatively  on  all  rays,  with  only  the  magnitude  and 
location  of  the  maximum  in  the  scattered  gas  contribution  changing  from  ray  to  ray. 
Figures  (3-5)  and  (3-6)  show  the  development  of  the  radial  profiles  with  increasing  0. 
Near  the  thrust  axis,  the  recovery  of  the  ambient  condition  is  slowest  because  the 
plume  is  densest.  The  density  overshoot  begins  at  about  0  =  15°,  and  increases  in 
magnitude  with  increasing  0,  until  peak  values  are  reached  for  values  of  0  sli  ltly 
larger  than  0  =  60°.  The  location  of  the  peak  moves  closer  to  the  nozzle  with  increas¬ 
ing  0.  For  larger  angles,  the  density  drops  rapidly  towards  the  ambient  value,  be¬ 
cause  the  exhaust  gas  density  is  so  low  that  no  significant  scattering  takes  place. 

The  results  predicted  by  the  present  theory  for  values  of  y  smaller  than  about 

_2  *  * 

10  ,  corresponding  to  r  =10,  must  be  treated  with  extreme  caution.  Since  r 

measures  the  radial  coordinate  in  units  of  exit  radii  (/?  is  about  0.3  in  the  present 

example),  phenomena  on  this  scale  may  be  strongly  influenced  by  the  details  of  the 

geometry  of  the  body  carrying  the  motor.  Moreover,  even  in  the  absence  of  a  finite 

body,  the  exhaust  gas  characteristics  in  the  vicinity  of  the  nozzle  are  not  described 
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accurately  by  Eq.  (3-2).  However,  the  solutions  for  larger  y  are  not  sensitive  to 

-2 

the  behavior  of  the  results  for  y  0(10  ). 

T  e  density  is  plotted  as  a  function  of  angle  at  several  radial  stations  in  Fig. 

(3-7).  The  ambient  conditions  are  recovered  progressively  later  in  the  expansion  as 
6  decreases,  since  the  region  near  the  thrust  axis  is  the  least  accessible  to  the  atmo¬ 
spheric  gases.  Figure  (3-8)  shows  the  corresponding  curves  for  the  windward  portion 
of  the  plane  containing  the  wind  and  thrust  vectors  in  the  non-aligned  case.  The  den¬ 
sity  distribution  is  qualitatively  similar  to  that  in  the  aligned  flow.  However,  the 
asymmetry  is  such  that  the  recovery  must  be  made  slower  on  the  leeward  side  (Fig. 
(3-9)),  because  of  the  fact  that  the  atmosphere  must  penetrate  through  the  entire 
windward  portion  of  the  plume  before  reaching  these  rays.  The  singular  ray  is  now 
at  9  =  90°,  on  the  leeward  side.  Although  the  model  predicts  some  scattering  along 
this  ray,  the  value  of  q^/W^  is  0(e  Aoc  /Kn),  which  is  less  than  10  in  the  present 
problem. 
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APPENDIX  A.  HYPERSONIC  LIMIT  -  HOMOGENEOUS  TERM 

The  homogeneous  term  in  Eq.  (3-7)  may  be  rewritten  (inverting  the  order  of 
integration)  as; 


2  r  t 


I ds'  h  I \ 8inJ  d?  JiF " 


-i\V2(  1-cos2  a  ) 

V  •>  \  / 


O  O 


Ho  ( W  cos  or ,  q  \  - 


(3-Al) 


A  ds  W  \X  +  V  s  j 
s' 


Here,  cos  a  is  given  by  Eq.  (3-6),  while  the  quantities  Tg.  a/,  W  ,  and  n/n^  are 
functions  of  x  +  vs' .  The  basic  assumption  is  analogous  to  that  employed  in  Section 
3.3,  that  the  dominant  contribution  to  the  angular  integrations  over  £  and  r?  comes  from 
the  neighborhood  of  cos  a  =  l,  because  of  the  large  factor  W2  in  the  exponential  in  the 
integrand  of  term  (3-Al).  By  letting  x  =  cos  £  be  the  first  integration  variable,  the 
exponential  factor  in  te’un  (3-Al)  becomes; 


-±w2  r-S- 


2  e  (r+;s-)2 


(3-A2) 


f(x)  =  1  -  1  x  cos  6  +•  v  1-  x  sinfl  cos  (p  -  <i>) 
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Although  the  quantity  multiplying  f(x)  in  (3-A2)  is  itself  a  function  of  x,  through  its 
dependence  on  the  variable  x  +  vs',  the  fact  that  it  is  large  and  positive  for  all  values 
of  v  means  that  the  contribution  to  (3-Al)  will  be  exponentially  small,  except  for 
values  of  v  corresponding  to  cos  or  =  1.  Examination  of  Eq.  (3-6)  shows  that  cos  o  =  1 
can  only  occur  for  v  =  r  .  Thus,  the  exponential  factor  (3-A2)  may  be  approximated  by:. 


1  w2  - L. - 

2  e  ,2 

(r  +  s' ) 


(3 -A3) 


W  is  a  function  of  r  (r  +  s')  in  Eq.  (3-A3).  The  method  of  steepest  descent  may  now 
0 

be  applied  to  (3-Al)  in  the  same  way  that  the  corresponding  integrals  for  the  inhomogeneous 

term  were  evaluated  in  Section  3.3.  The  result  implies  that  the  remaining  terms  in 

the  integral  are  to  be  evaluated  at  x  =  cos  9,  tj  -  <t>.  When  these  values  are  inserted 

into  the  expression  (3-6),  then  cos  0=1,  as  required  for  consistency,  provided  that 

s'  +  r  is  positive.  Where  s'  +  r  is  negative,  cos  o  =  -1,  and  the  contribution  to 

(3-Al)  is  again  exponentially  small  in  W  ,  this  time  because  of  the  fact  that  the  func- 

-W^  /2  e 

tion  H„(-W  ,q  )  is  0(e  ),  if  now.  the  substitution  r'=  s'+  r  is  made  in  the  integra- 

o  0  0 

tion  over  s'  which  remains  in  Eq.  (3-Al),  then  the  homogeneous  term  in  Eq.  (3-11)  is 
recovered. 
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APPENDIX  B.  PARABOLIC  INTEGRATION  FORMULAE 
WITH  VARIABLE  INTEGRATION  STEP 


Consider  the  following  mesh  of  three  points  x 2,  and  x^: 


xl>  x2  =  X1  +  hl  ’  x3  =  X2  +  h2 


and  the  following  integrals  of  a  function  f(x): 


f(x)  dx 


(3-B1) 


*2,3  =  J  rWdX 


(3-B2) 


f(x)  dx 


(3-B3) 


Obviously 


^13  ^1  2  ^  ^2  3 


(3-B4) 
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The  objective  here  is  to  derive  approximate  expressions  for  these  integrals,  by 
use  for  f(x)  of  the  second-order,  lagrange  interpolating  polynomial:. 


f(x)  = 


(x  -  x2)  (x  -  x3)  -  iqr2  (X  -  Xl)  (X  -  X3) 

+  MVV)  (x'xi)(x'x2) 


(3-B5) 


where  f ^  =  ffx^),  and  so  on. 


jlK-x2)(x-x3)dx=i(x»-^)-i(x2+x3)(x2-xf;;+x2x3lx2-Xi) 


6  (x2'xlJ2  (3x3-x2-2xl) 


=  f(3h2+2hl)  ' 


(3-B6) 


By  interchanging  x^  and  x9  in  Eq.  (3-B6),  one  obtains: 


j  (x-xl)(x-x3)dx'l(xl-x2)2  (3X3-X1'2X2) 


=  -  -r  <3h2  *  hji 


(3-B7) 
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Finally,  by  adding  Ecs.  (3-B9)  and  (3-B10),  one  obtains: 


h  +  h  2h  -  h 
_£ - i  — ± - £  f  + 

6  hj  1 


(h2  +  hl) 

“7*T 


2h  -  h 

f  +  - - - i  f 

*2 


(3-Bll) 


The  error  involved  in  the  approximate  integration  formulae  (3-B9)  to  (3-Bll)  is 

4 

of  the  order  of  h  .  For  hj  =  h2>  formulae  (3-B9)  and  (3-B10)  reduce  to  the  well-known 
formula  of  Adams-Moulton,  while  formula  (3-Bll)  reduces  to  Simpson's  rule. 
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4.  CLASSICAL  CALCULATIONS  OF  ROTATIONAL  AND  VIBRATIONAL 
EXCITATION  IN  ENERGETIC  ATOM  -DIATOMIC  HYDRIDE  COLLISIONS 


4.1  Introduction 

Since  1960,  literally  hundreds  of  calculations  of  inelastic  energy  transfer 
in  simple  atom -diatom  or  diatom -diatom  collisions  have  been  published.  The 
overwhelming  bulk  of  this  work  has  been  pe  rformed  while  using  either  the  rigid 
rotor  approximation  to  study  only  translator ial  -rotational  energy  transfer  or  a 
collinear  collision  scheme  to  isolate  translational -vibrational  energy  transfer. 

Recently,  a  number  of  workers  have  published  studies  that  treat  the 

problem  of  coupled  vibrational  and  rotational  energy  transfer  in  simple  collisions. 

These  efforts  include  the  classical,  two-dimensional  models  of  Benson  and 

Berend,  Kelley  and  Wolfsberg,  and  Bergeron  and  Chapuisat,  ^^which 

aim  primarily  at  elucidating  the  effects  of  rotational  inelasticity  on  collisionally 

induced  vibrational  excitation  or  de -excitation.  Semiclassical,  three-dimensional 

treatments  of  this  same  problem  have  been  presented  recently  by  Wartell  and 
(20)  (?U  (22) 

Cross,  '  Stallcap  and  Sorensen;  '  Other  studies  of  interest  incluue  the 

quantum  mechanical  treatment  of  the  He-H„  collision  system  by  Eastes  and 

(23)  ^ 

Secrest, '  ’  which  uses  the  He-H„  interaction  potential  calculated  by  Gordon  and 

(24)  ^  (25) 

Secrest, '  '  and  the  semiclassical  treatment  of  the  Ar-H2  system  by  Shin, '  ' 

which  uses  an  interaction  potential  based  on  molecular  beam  data.  Doll  and 

Miller^26^  have  also  treated  the  He-H2  problem  by  using  classical  S-matrix 

(27) 

techniques.  In  addition,  Thompson'  has  been  able  to  calculate  reasonable 
rotational  and  vibrational  energy  transfers  for  a  number  of  atom -diatom  systems 
by  using  Monte  Carlo  classical  trajectory  techniques  and  semiempirical  valence - 
bond  interaction  surfaces. 
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The  focus  of  the  present  work  is  the  effect  of  vibrational  inelasticity 
on  the  rotational  excitation  of  such  diatomic  hydrides  as  HF,  HC1,  and  OH 
in  energetic  collisions  with  atomic  species.  Classical,  three-dimensional 
calculations  of  rotational  excitation  probabilities  for  both  rigid  and  nonrigid 
diatomics  in  collisions  with  several  atomic  species  are  presented  here,  and 
the  effect  of  the  rigid  rotor  approximation  is  evaluated. 

4.2  Interaction  Potentials 

The  collisions  of  interest  have  center -of -mass  energies  in  the  range  of 
1  to  5  eV.  Such  interactions  occur  well  up  on  the  repulsive  core  of  reasonable 
intermolecular  potentials  and  are  insensitive  to  the  weaker  "long-range"  portions 
of  the  potentials.  This  insensitivity  is  exploited  by  adopting  hard-sphere  poten¬ 
tials  for  the  colliding  species.  Each  diatomic  is  modeled  with  two  interpenetrating 
hard  spheres,  while  monatomic  species  are  represented  by  a  single  sphere.  The 
diatomics  are  either  held  rigid  along  their  line -of -centers  at  a  distance  between 
sphere  centers  fixed  by  the  equilibrium  ground-state  internuclear  distance  of  the 
hydride  in  question  or  are  attached  along  the  line-of-centers  by  a  one -dimensional 
harmonic  oscillator. 

Sphere  sizes  for  bound  and  unbound  atoms  are  estimated  from  the  high- 

energy  molecular  beam -scattering  data  of  the  Amdur -Jordan  group  at  MIT  or 

of  the  Leonas  group  at  the  Moscow  State  University.  The  sphere  sizes  selected 

represent  the  effective  range  of  one  eV  potential  for  the  species  in  question  and 

/og\ 

are  determined  by  methods  outlined  in  earlier  work  on  atom-H90  collisions. ' 

Equilibrium  bond  distances  for  the  hydride  molecules  are  taken  from 
Herzberg,^^  and  are  shown  in  Table  4-1, along  with  hard-sphere  radii  for  mole¬ 
cular  atoms  and  monatomic -collision  partners.  Sphere  radii  for  O,  H,  Ar  and 

He  are  taken  from  Ref.  (28).  The  F-atom  radius  is  based  on  data  of  Belyaev, 

(30)  9  30 

et  al. ,  '  which  determined  a  potential  of  V(Ar-F)  of  506/r  *  ,  where  \  is 

expressed  in  eV  and  r  in  A.  The  Cl  radius  is  estimated,  since  no  pertinent 
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molecular  beam  data  could  be  located.  The  same  radius  and  internuclear 

35  37 

separation  is  used  for  cases  involving  the  isotopes  Cl  and  Cl.  All  other 
atoms  are  assumed  to  have  only  the  mass  of  their  predominate  isotope. 


Table  4-1.  Molecular  Parameters 


Molecular  Atom 

Hard-Sphere  Radius  (A) 

H 

0.65 

0 

0.88 

F 

0.83 

Cl 

1.18 

Monatomic  Collision 

Partner 

Hard -Sphere  Radius  (A) 

He 

0.63 

O 

0.85 

Ar 

1.12 

Equilibrium  Bond 

Average  Geometric 

Diatomic 

Distance  =  (A) 

O 

Cross  Section  (A  ) 

OH 

0.9706 

3.06 

HF 

0.9171 

2.80 

HC1 

1.2746 

4.92 

4.3  Collision  Mechanics 

The  choice  of  a  hard  sphere  interaction  potential  reduces  the  analysis  of 
the  collision  process  to  a  study  of  the  configuration  at  the  moment  of  in  pact.  The 
calculation  is  most  conveniently  performed  in  the  precollision  rest  frame  of  the 
molecule.  The  rotational  and  vibrational  degrees  of  freedom  of  the  molecule  are 
assumed  to  be  initially  unexcited.  Thus,  the  configuration  is  fixed  in  an  inertial 
frame  at  the  collision  instant. 


Now  consider  the  impact  of  the  collision  partners  q  with  the  molecular 
atoms  denoted  by  the  subscripts  1  and  2  (Fig.  4-1).  The  conservation  of  momen 
turn,  angular  momentum,  and  energy  are  given  respectively  by 


m_.  v_  =  m_  v'  +  m.  y’ 


q-q 


q  -q 


1  *1 


m„  v 


2  -2 


(4-1) 


m  r  xv  =  m  r  x  v*  +  m,  r,  x  v'  +  m_  r„  x  v' 
q-q  -q  q  -q  q  1-1  -1  2-2  -2 


(4-2) 


1  2  1  ,2 

v  =  rm  v' 

2  q  q  2  q  q 


12  12 
=  2miV'l  +  2  “2  ^2 


(4-3) 


Figure 4. 1  -  Coordinate  systems  used  to  describe  collision. 
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The  primed  quantities  denote  post  collision  values  of  the  atom  velocities 
va(Q!  =  l,2,q),  while  r^  denotes  the  position  of  the  center-of-mass  of  the 
«  atom  with  respect  to  the  center-of-mass  of  the  molecule  at  the  instant  of 
collision. 

In  order  to  proceed  further,  it  is  necessary  to  consider  the  relationship 
between  the  hard-sphere  model  and  a  continuous  molecular  interaction  poten¬ 
tial  <f>(r)  °f  the  form 


*<r>  *  *1  (rlq)  +  *2  (r2q) 


2  M12 


2  2 
w  $ 


(4-4) 


r 


aq 


a  =  1,2 


(4-5) 


Here  p12  is  the  reduced  mass  of  the  molecular  atoms,  w  is  the  frequency  oi 
vibration,  and  £  is  the  amountby  which  the  distance  between  the  molecular 
atoms  has  departed  from  its  equilibrium  value  (the  "spring"  deflection).  For 
collisions  sufficiently  energetic  to  penetrate  well  up  on  the  repulsive  core  of 
either  atomic  potential,  ^Q,(raq)>  the  local  length,  f .  over  which  the  interaction 
occurs  may  be  characterized  by 


lr“1=r»)  •  <4'6) 

The  radius,  r  ,  is  chosen  so  that  (r  )  is  in  the  interaction  energy  range  under 

°  28j  a  0 

consideration.  '  '  The  interaction  time,  t,  is  approximately 


T 


n 


V 

^4 


(4-7) 
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The  hard-sphere  model  may  be  regarded  as  the  limiting  form  of  <f>( r)  as 
f/r  approaches  zero,  provided  that  two  further  conditions  are  satisfied. 
First,  the  rotational  excitation  must  be  such  that  |  t  «  1.  The  angular 
velocity  Q  is  defined  so  that  the  impulsively  generated  angular  momentum  of 
the  molecule  about  its  center  of  mass  is  in,  where  I  is  the  moment  of  inertia. 
Under  these  conditions,  the  magnitude  of  the  impulse  W  transmitted  to  the 
collision  partner  q  is  directed  along  the  line  connecting  the  centers  of  the  two 
colliding  atoms.  Thus, 


=  n  W 


(4-8) 


where  is  the  unit  normal  to  the  molecular  atom. 

The  manner  in  which  the  impulse, is  divided  between  the  molecular 
atoms  is  determined  from  the  final  condition  that  must  be  satisfied  for  the 
hard-sihere  model  to  be  employed.  If  the  frequency  of  the  vibration  is  suffi¬ 
ciently  -ow  so  that  lot  «  1,  then  the  impulse  is  transmitted  before  the  "spring" 
(31) 

can  deform.  Since  the  impulse  must  be  transmitted  to  the  second  molecular 
atom  via  the  vibrational  mode,  the  instantaneous  post-collision  state  of  the 
molecule  n  ust  be  such  that  the  molecular  atom  experiencing  the  impact  absorbs 
the  component  of  the  impulse  directed  along  the  spring  axis,  k.  (See  Fig.  4-1.) 


-n,k'\V  =  m,v,  •  k 
-1"  1-1  - 


(4-9) 


If,  on  the  other  hand,  the  vibrational  frequency  is  sufficiently  high  so  that 
lot  »  1,  then  the  molecule  will  behave  like  a  rigid  rotor,  and  Eq.  (4-3)  must 

(3D 

be  replaced  by  the  condition 

£<t)  =  0  .  (4-10) 
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When  cur  is  neither  large  nor  small,  the  use  of  the  hard-sphere  model  cannot 
be  justified  on  the  basis  of  classical  mechanics.  However,  in  order  to  apply 
the  classical  approximation  to  a  molecule  containing  a  vibrational  degree  of 
freedom,  the  classically  computed  vibrational  excitation  must  be  much  larger 
than  the  energy  hcu  needed  to  excite  that  degree  of  freedom.  When  this  is  not 
the  case,  the  rigid  rotor  assumption  [  Eq.  (4-4)]  may  still  provide  a  useful 
approximation  to  the  rotational  excitation,  even  though  the  condition  cut  «  ] 
is  not  satisfied.  This  point  is  discussed  further  in  Sec.  4.7. 

The  system  of  Eqs.  (4-1)  through  (4-4)  is  now  used  to  determine  the  excita¬ 
tion  of  the  internal  degrees  of  freedom  in  a  single  collision  for  both  the  rigid 
rotor  (cut  »  1)  and  the  harmonic  oscillator  model  (cur  «  1).  To  this  end,  it  is 
convenient  to  replace  the  molecular  velocities,  v*  ,  by  a  center-of-mass 
velocity,  vc^,  an  angular  velocity,  ft,  and  a  spring  deflection  velocity,  w*k. 
These  quantities  are  defined  by  the  conditions 


(,ml  +  m2)  *CM  =  ml*l  +  v2“2  ’ 


Ift  =  m,a,k  xv',-  m-a^k  x  v'_  . 


1*1“  *2)  **  =W 


(4-11) 

(4-12) 

(4-13) 


thus, 


m2w 

v',  =  V-,..  +  ft  x  a,k  +  - ; - 

-I  -CM  1-  m  +  m2 


k  , 


(4-14) 


miw 

vk  =  -  ft  x  a~k  -  — — t - k  . 

-2  -CM  -  2-  +  m2 


(4-15) 
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By  substituting  Eqs.  (4-14)  and  (4-15)  into  Eqs.  (4-1)  through  (4-3),  the  conserva¬ 
tion  laws  take  the  form 


( m,  +  m0\  v„.,  =  m  (v  -  v'  \ 
VI  2)  -CM  q  \-q  ~qj 

(4-16) 

mr  xiv  -  V  \  =Ifi  , 

(4-17) 

q-q  v'q  -qj 

1  2  1  2'  1 

r-m  v  =  Trm  v  +  r/m,+  m,| 

2  q  q  2  q  q  2  V  1  2 

2  L  lrr>2  ^  1  2* 

i  VCM  2If2  2m12W 

(4-18) 

Equations  (4-2)  and  (4-6).  together  with  the  expression 

Hlw*  =  mi XCM  *  k+  M12w*  (4-19) 

for  the  oscillator  model,  and  the  condition 

w*  =  0  (4-20) 


for  the  rigid  rotor  model,  constitute  the  formulation  of  the  problem. 

The  above  system  of  equations  is  readily  solved  by  using  Eq.  (4-8)  to 
eliminate  from  Eqs.  (4-16)  and  (4-19).  Thus.  v^^j,  2-  an(^  w*  can  be  ex¬ 
pressed  as  functions  of  the  impulse  W  as 


n,  W 
-1 


-CM  mi+  m2 


/  r  x  n. )  W 

v'q  -ij 


n  = 


* 

W  = 


(Si  •  &)  w 


m. 


(4-21) 


(4-22) 


(4-23) 
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The  last  equation  is  replaced  by  w*  =  0  for  the  rigid  rotor  model.  Substitution 
of  Eq.  (4-22)  into  the  energy  conservation  equation  yields  the  following  results 
for  W  by  using  the  harmonic  oscillator  model 


W  =  - 


2v  •  n. 


n,  •  k 


\  2 


-1  |  ■•1  ~  \ 

M  +  “12 


/  r  x  n,  \ 

l~q  -1; 


(4-24) 


where 


r  =  a.k  +  /R.  +  R  )  n. 
-q  1-  \  1  qj  1 


and 


m  ■m 
D  a 

n  = - ■ — ,  mn  =  m,  +  m0  . 

qD  mn  +  m  D  1  2 

M  D  q 

Here,  R  and  R  are,  respectively,  the  hard-sphere  radii  for  the  molecular 
t  q 

atom  experiencing  the  collision  and  the  collision  partner.  The  corresponding 
expression  for  the  rigid-rotor  model  is 


W  =  - 


2v  •  n, 
~q  -*-1 


_x  irqXnl) 

(MqD)  +  I 


21 


(4-25) 


Expressions  (4-21)  through  (4-25)  constitute  the  required  solution.  Ignoring 
those  incident  configurations  that  lead  to  multiple  impacts  (e.g. ,  certain  col- 
linear  collisions)  as  statistically  unimportant  compared  with  those  that  do  not, 
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the  rotational  and  vibrational  excitation  energies  experienced  in  an  individual 
collision  are  given  by 

EV=2M12W*2  ’  (4-26) 

En  =  |l«2  .  (4-27) 

These  expressions  follow  immec'r  telv  from  the  assumption  that  the  deformation 
|  is  small  compared  with  a1  +  a 2,  s  ,  that  the  post  collision  oscillations  can  be 
expressed  in  terms  of  a  Lagrangian,  Ly>  of  the  form 

Ly  =  2M12(W*)2  “  2M12w2^2  •  (4-28) 

The  modification  of  these  results  because  of  the  presence  of  a  steep,  con¬ 
tinuous  potential  of  the  type  considered  above  could  be  calculated,  in  principle, 
using  multi-time  scale  perturbation  techniques.  For  example,  such  a  calcula¬ 
tion  could  be  employed  to  determine  the  effect  of  a  slow  rotation  or  vibration  of 
the  molecule  about  the  center  of  mass  on  the  trajectory  of  the  collision  partner 
as  it  passed  through  the  potential,  of  one  of  the  molecular  atoms. 

The  rotational  and  vibrational  excitation  would  be  altered  by  an  amount  of  the 
order  of  |Sjr  by  this  phenomenon,  which  would  permit  vibrational  excitation  of 
this  magnitude  to  occur  even  tor  collisions  perpendicular  to  the  molecular  axis. 

Similar  procedures  could  be  employed  when  the  vibrational  period  is  much 
smaller  than  the  interaction  time.  However,  it  is  not  clear  that  such  elaborate, 
classical  calculations  of  vibr  itional  effects  should  be  performed  until  it  is 
evident  that  a  classical  description  of  this  channel  is  appropriate. 


1-64 


4.4  Kinematics  and  Probability  Model 

The  colliding  atom  is  assumed  to  be  drawn  from  a  uniform  beam  of  such 

atoms  with  common  velocity  vector  v  .  The  diatomic  molecule  is  assumed  to 

~q 

have  a  random  orientation  with  respect  to  its  center-of-mass  at  the  time  of 
impact.  In  what  follows,  the  (i,  j,  k)  coordinate  system  is  fixed  in  the  molecule 
with  origin  at  the  center-of-mass  and  the  axis-of-symmetry  in  the  z-direction 
(see  Fig.  4-1),  so  that  the  beam  is  viewed  as  randomly  oriented.  By  sym¬ 
metry,  Vq  can  be  assumed  to  lie  in  the  y-z  plane  for  all  collisions,  making  an 
angle  0  with  the  z  axis,  with  differential  probability  density,  sin  0d0.  The 
location  of  the  collision  partner's  center  prior  to  impact  is  given  by  the  impact 
parameter  b  and  an  angle  ft,  which  are  polar  coordinates  in  a  plane  orthogonal 
to  v  .  Since  this  distribution  is  uniform,  the  differential  probability  density  is 
b  db  d/3.  The  third  location  parameter,  D,  the  initial  distance  apart,  is,  of 
course,  arbitrary.  A  second  coordinate  system  (i*,  j*.  k*),  centered  in  the 
colliding  atom,  is  defined  by  rotation 

i*  '=  i.  j*  =  j  cos  0  -  k  sin  0,  k*  =  j  sin  0  +  k  cos  0  ,  (4-29) 

so  that  the  beam  axis  is  k*.  It  follows  that  the  initial  location  of  the  atomic 
center  is 

p  =  Dk*  +  b  (cos  ft  if  +  sin  /3  j*)  . 

=  (D  sin  0  +  b  sin  ft  ccs  0)  j  +  (b  cos  /3)  i  , 

+  (D  cos  0  -  b  sin  ft  sin  0)  k  .  (4-30) 

The  point  of  impact  with  one  of  the  molecular  atoms  (say  atom  1)  is  deter¬ 
mined  next.  Let  the  vector  between  atomic  centers  initially  be  c  -  p  -  a,k. 

At  the  time  of  impact,  the  colliding  atom  has  traveled  a  distance  p  in  the 
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direction  -k*.  Let  the  final  vector  between  centers  be  (R,  +  R  )n'  ,  where 

|  q  «*** 

n'  is  a  unit  vector  relative  to  atom  1.  It  follows  that 

(Rl+Rq)^  "S.-PS*  •  <4'31> 

which  defines  both  p  and  n'  .  The  local  angles  9'  and  a'  are  found  from 
cos  9’  =  n'  •  k  and  cos  a'  =  n'  •  k*.  The  solution  to  these  equations  yields 


sin  a'  =  J  ^  sin  6  +  b  sin  0j2  +  (b  cos  +  R^j  ,  (4-32) 

p  =  D  -  a^  cos  0  -  ^Rj  +  R^j  cos  a'  ,  (4-33) 

cos  9'  =  cos  9  cos  a'  -  sin  9  ^  sin  0  +  b  sin  (3\/(R^  +  R  \  .  (4-34) 


where,  if  the  r.h.s.  of  the  first  equation  is  greater  than  unity,  there  is  no 
collision.  In  the  case  of  impact  with  atom  2,  one  must  replace  Rj  by  R2,  and 

al  by  ~a2" 

Finally,  let  p  ,  p2  denote  the  distance  p  from  impact  with  atom  1  and 
atom  2.  respectively.  Then,  the  actual  impact  will  occur  on  atom  1  if  p  <  p^, 
and  on  atom  2  otherwise  (which  choice  is  independent  of  Dl.  Also,  given  the 
angles  0  and  p,  a  maximum  feasible  impact  parameter,  b  .  is  obtained  by 
setting  sin  a'  =  1. 


4 . 5  Calculation  of  Rotational  and  Vibrational  Energy  Distributions 

For  each  collision,  the  rotational  and  vibrational  energies  are  given  by 
Eqs.  (4-26)  and  (4-27),  where  for  the  rigid  rotor,  Ey  =  0.  In  the  following 
treatment,  these  energies  are  normalized  by  the  center-of-mass  energy, 
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ECM  ~  -1/2  MqD  Vq-  Assuming  for  definiteness  that  the  impact  occurs  on  atom  1, 
in  the  vibrational  case  the  formulas  simplify  to 


ER'/ECM  =  E1  sin2  e'  cos"  a'  ’  Ev//ECM  =  E1  cos2  cos2  a' 


(4-35) 


Thus,  the  maximum  energies,  and  E^.,  occur  for  normal  collisions  at 
O'  -  tt/2  and  O'  =  0,  respectively.  The  maximum  normalized  value  in  either 


case  is 


E  =  4 
1  2 

mi 


MqDA'l2 


(4-36) 


which  is  a  function  of  the  masses  only  and  never  exceeds  unity.  It  is  also 
easily  seen  that  E^  is  larger  than  E 2  (corresponding  to  impact  on  the  heavier 
atom  2). 

In  the  rigid  rotor  case,  the  formula  (for  atom  1)  becomes 


^CM  ’  4 


^qD  ^12 


1  +  —  sin2  O'  )  sin  9'  cos2  o”  . 

mi  / 

1  '  (4-37) 


In  this  case,  the  maximum  normalized  energy  is  a  function  of  O' ,  as  well  as  the 
masses,  and  attains  the  value  unity  for  a  normal  collision  with 


2  m  m  +  m 

sin  O'  =  —  1  +  — - - 

mq  \  m2 


(4-38) 


(if,  and  only  if,  the  r.h.s.  is  less  than  unity). 
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2 

This  condition  is  equivalent  to  >1,  and,  if  violated,  the  maximum  value 

is  identical  to  the  vibrational  case  (namely,  at  0'  =  7r/2).  Again,  impact  on 
the  heavier  atom  clearly  leads  to  higher  energies. 

The  probability  distributions  for  the  normalized  energies  are  computed  by 
a  Monte  Cav’o  technique  for  each  of  the  collision  pairs  in  Table  4-1,  and  for 
both  the  rigid  rotor  and  vibrational  models.  Pseudorandom  numbers  are  used 
to  generate  sample  collisions  with  the  collision  parameters  (0,  (3,  b)  drawn 
from  a  probability  density  with  differential 

b  sin  0  dO  dbd/3  .  (4-39) 

9 

This  amounts  to  generating  uniform  random  variables  (cos  0,  (3,  b  )  and  deter¬ 
mining  the  local  collision  angles  (O' ,  cr)  or  (0"  ,  or  )  by  the  formulas  of  the 
previous  section.  By  considering  symmetry  and  realizability,  the  collision 
parameters  are  restricted  to  the  intervals  0  <0  <n,  6  </3  <w,  and 
O<b<bmax(0,  p). 

From  the  sample  collisions,  a  frequency  histogram  of  E/E^M  is  con- 

r  v  v 

structed  for  each  of  the  three  energies,  E^,  E^,  and  Ey.  This  histogram  is 
converted  to  an  approximate  probability  density  P(x)  by  the  definition 

P(x)  =  prob  -|x.  <  E/Ecm  <  x.  +  Ax  j-  /Ax  ,  (4-40) 

for  x.  <  x  <  Xj  +  Ax.  In  all  cases,  20  subintervals  are  used  on  the  interval 
0  <  E/ECM  -  (E/'ECM)max  ,  so  that  Ax  -  0.05  E^  for  the  vibrational  model,  and 
Ax  =  0.  05  for  the  rigid  model.  The  total  number  of  sample  points  used  is 
50,000.  The  results  are  displayed  in  Figs.  4-2  through  4-10. 
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Figures  4. 2  through  4. 10  -  Normalized  probability  distributions  for  the  inelastic 

excitation  of  the  diatomic  hydrides  OH,  HF.  and  HC1 
in  collisions  with  He.  O,  and  Ar,  Solid  histogram 
represents  E^/ECM,  and  drshed-dotted  histogram, 
Er/ECM’  for  the  harmonic  oscillator-atom  colli¬ 
sions.  Dashed  histogram  represents  E^/E^j  f°r 
rigid-rotor-atom  collisions.  Specific  sets  of  colli¬ 
sion  partners  are  indicated  for  each  figure. 
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Figure  4.2  -  O  +  HF 
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Figure  4.4  -  He  +  HF. 


Figure  4.5-0+  HC1. 
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Figure  4.6  -  Ar  +  HC1 
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Figure  4.8  -  O  +  OH. 
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35  37 

Separate  runs  were  made  for  atomic  collisions  with  H  Cl  and  H  Cl. 
However,  the  probability  distributions  showed  virtually  no  isotope  effect,  and 
only  the  results  for  H^5Cl  are  displayed.  Table  4-II  gives  the  sample  means 
in  the  following  format.  The  top  number  for  each  collision  pair  represents  the 


Table  4-n.  Normalized  Mean  Excitations 


Diatomic 

Atom 

He 

O 

Ar 

er/ecm 

0.127 

0.112 

0.113 

^CM 

0.083 

0.069 

0.072 

OH 

</ecm 

0.036 

0.029 

0.029 

pl 

0.283 

0.296 

0.309 

Er/ecm 

0.127 

0.106 

0.103 

HF 

Er/Ecm 

0.083 

0.063 

0.062 

Vecm 

0.035 

0.026 

0.025 

pl 

0.300 

0.319 

0.329 

Er/ECM 

0.077 

0.059 

0.055 

eXm 

0.044 

0.  029 

0.030 

HC1 

ev/ecm 

0.023 

C.  015 

0.014 

pl 

0.186 

0. 195 

0.210 
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mean  fraction  of  the  center-of-mass  collision  energy  which  the  rigid  rotor-atom 
collision  imparts  to  the  diatomic’s  rotational  modes  (E^/E^M).  The  second 
number  for  each  pair  is  the  mean  fraction  imparted  to  the  rotational  modes  in 
the  harmonic  oscillator-atom  collision  (E^/E^^).  The  third  number  repre¬ 
sents  the  mean  fractional  energy  in  the  vibrational  mode  for  the  harmonic 
oscillator-atom  collisions  (Ey/Ec^).  The  sample  probability  px  of  impacting 
on  atom  1  (the  H  atom)  rather  than  atom  2  is  also  shown  as  the  fourth  number. 
Estimates  of  the  maximum  relative  error  in  the  plotted  values  of  P(x)  can  be 
made  by  using  a  normal  approximation  to  the  multinomial  distribution  for  a 
histogram.  These  estimates  typically  range  from  2%  for  P(0)  to  100%  for 

P(E  ).  The  relative  errors  in  the  means  and  p1  are  all  within  2%. 

'  max'  1 

4.  6  An  Approximate  Distribution  for  the  Vibrational  Model 

The  simplicity  of  the  energy  expressions  in  the  vibrational  case  leads 
naturally  to  an  approximate  analytic  expression  for  P(x)  which  is  quite  accurate . 
First,  assume  that  the  probabilities  p^.  p9  of  impacting  atom  1  or  atom  2. 
respectively,  have  been  accurately  estimated  by  Monte  Carlo  methods,  and  then 
restrict  attention  to,,  say,  atom  1.  The  probability  distributions  for  the  random 
variables  (O'  ,  of  )  are  derived  from  those  for  (0,  b,  /?),  and  then  transformed 
analytically  to  ohtain  the  distributions  for  the  energies.  However,  to  make  the 
transformation  tractable,  a  certain  geometric  approximation  is  introduced. 

From  physical  considerations,  one  expects  the  distribution  of  hits  on  the 
surface  of  atom  1  to  be  random;  that  is,  cos  0'  is  uniformly  distributed..  This 
conjecture  is  verified  both  by  simulation  and  manipulation  of  the  equations  of 

9 

transformation.  Similarly,  inspection  of  the  equation  for  sin"  a'  leads  to  the 

2 

conclusion  that  it  (and  hence  cos  Q"  )  is  also  uniformly  distributed  indepen¬ 
dently  of  O'  . 
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Both  these  results  ignore  the  effects  of  shielding  by  atom  2  on  some 
trajectories  that  would  otherwise  impact  atom  1.  Such  shielding  is  of  two 
types  (see  Fig.  4-11).  In  the  first  type  (total  shielding),  a  range  of  angles 
0^<  6'  <  rr  is  totally  inaccessible.  In  the  second  type  (partial  shielding),  for  a 
range  of  angles  0'^,  <  0'  <  0'^  the  probability  of  impact  is  reduced,  and  the  dis¬ 
tribution  of  a'  depends  on  0' .  For  0  <  0'  <  0'T  atom  2  has  no  effect.  The 
approximation  made  is  to  ignore  partial  shielding,  while  taking  account  of  total 
shielding.  The  formula  for  0'^  is  easily  obtained  from  the  geometry, 


Figure  4.11-  Definition  of  the  angles  0'^  and  0'^ . 


1-81 


Now,  still  restricting  attention  to  atom  1,  let 


o 

x.^  =  cos  a'  ,  x2  =  cos  6’  .  (4-43) 

Then, 

2  ( .  2  N 

yi  X1X2  ’  y2  X1V  "X2  )  • 

and  (Xj,  x2)  are  uniform  random  variables  with  ranges  0  <  Xj  <  1  and 
-dj  =•  cos  0'K  <  x2  <  1.  (Note:.  0^  <  tt/2  for  the  collision  pairs  considered.) 

It  is  now  straightforward  to  derive  the  distribution  of  (y^, y2).  Using  the 
Jacobian  of  the  transformation  J  =  the  joint  distribution  is  found  to  be 

i  i  r  i-1/2 

dxtdx2  =  yi  (yi +  y2) j  dyidy2  •  <4-44> 

from  which  the  marginal  distributions  for  y^  and  y9  are  obtained  by  integration. 
Denote  the  joint  density  function  by  f(y^,  y2).  In  order  to  do  the  integrations,  it 
is  necessary  to  specify  the  region  of  admissible  (y^yglT  outside  of  which 
f(yry2)  sO]  .  Since  the  Jacobian,  J,  equals  zero  along  the  line  x2  =  0,  the 
integrals  must  be  calculated  separately  over  different  regions,  according  to  the 
sign  of  x2<  By  using  the  bounds  on  (Xj.Xg),  these  regions  are  and  as 
indicated  in  Fig.  4-12;  i.e.. 


(  (  f(iyry2)ds'idy2  =  1  • 

(4-45) 

jJ  \  / 
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2  ^  f tyl’ y2jdyl  “  1^7 00S  h_1  (y2l/2)  • 

o 


f(yry2)dyi 


(4-47) 


V. 


if  y2  <  i  -  <5j 


The  means  are  found  to  be  =  l/6  [  (1  +  <5^  )/(l  +  )]  ,  and  y2  =  1/2  -  yy 

where  the  last  relation  follows  from  the  identity  y^  +  y2  =  x^. 

Now  consider  impacts  on  atom  2.  Exactly  the  same  reasoning  leads  to 
identically  the  same  distributions  if  is  replaced  by  62  =  cos 
(Note:  0^  <  ir/2.)  The  distributions  on  the  two  atoms  separately  can  be  folded 
together  with  the  empirical  probabilities  p.^  and  p2  to  yield  the  total  density 


E 


V 


lE 


CM 


v 

) 


plfl 


ev/ecm 


|/E1> 


//E1+ P2fl  (^EV/Ecm)^E2’62  |/E2 


(4-48) 


and  similarly  for  E^.  replacing  f  by  f2-  (Recall  that  dyj  =  d(Ey/ECM)/E1< ) 

In  Fig.  4-13,  this  approximate  density  is  compared  with  the  Monte  Carlo 
results  for  one  case  and  is  observed  to  be  within  the  statistical  error.  The 
approximate  and  sample  means  differ  by  less  than  2%.  (Note:  A  similar 
approximation  can  be  constructed  for  the  rigid  rotor  model,  but  would  involve 
numerical  integration  for  the  densities.) 
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Figure  4.13-  Comparison  between  approximate  and  empirical  distributions  for 
O  +  HF^within  the  /ibrational  model.  Dashed  histogram  repre¬ 
sents  Ey/E^j^j,  and  dotted  histogram,  Er/Eqj^,  and  smooth 
curves  represent  approximate  solutions  for  these  distributions 


4.7  Differential  Energy  Transfer  Cross  Sections 


The  probability  distributions,  P(E/Ecm),  shown  in  Figs  4-2  through  4-11, 
can  be  converted  to  differential  energy  transfer  cross  sections  describing  the 
harmonic  oscillator  rotational  or  vibrational  excitation  or  the  rigid  rotor  rota¬ 
tional  excitation  as  a  function  of  incident  center-of-mass  kinetic  energy  by  the 
relation, 


(e/ecm) 

9E 


ghsP(E/ECM) 

ecm 


(4-49) 


where  is  the  geometric  hard-sphere  cross  section  for  the  appropriate  atom- 
diatom  pair,  and  E  =  E^,  E^,  or  E^,  respectively. 

Approximate  hard-sphere  cross  sections  for  the  system  studied  can  be 

calculated  by  the  following  procedure:  The  average  geometric  cross  sections 

for  the  hard-sphere  diatomic  models  are  listed  in  Table  4-1.  These  cross- 

sectional  areas  correspond  to  spherical  areas  of  projection  with  radii  equal  to 

0.99A,  0.94A,  and  1.25AforOH,  HF,  and  HC1,  respectively.  Normal  hard- 

2 

sphere  collision  cross  sections  (in  units  of  A  )  for  each  atom-diatom  pair  can 
then  be  constructed  from  the  relation 


ahs 


(4-50) 


where  d  ^  is  equal  to  the  effective  spherical  radius  for  the  appropriate  diatomic 
as  listed  above,  plus  the  hard-sphere  radius,  R^,  of  the  appropriate  atomic 
collision  partner,  as  listed  in  Table  4-1. 
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4.8  Discussion 


The  average  energy  transfer  results  listed  in  Table  4-II,  and  the  energy 
transfer  probability  distributions  shown  in  Figs.4-2  through  4-10  can  be  inter¬ 
preted  by  noting  that  energy  transfer  in  the  classical,  hard-sphere  collision 
model  used  in  this  study  depends  only  on  the  reduced  mass  of  the  atom-diatomic 
hydride  system  and  on  geometrical  constraints  that  are,  in  turn,  dependent  on 
hard-sphere  radii  and  bond  distances.  These  geometrical  constraints  arise 
partly  from  the  fact  that  some  impact  points  on  the  diatomic  molecule  are 
inaccessible  from  certain  directions  due  to  shielding  by  the  second  molecular 
atom,  while  other  impact  points  are  completely  inaccessible  due  to  the  finite 
size  of  the  colliding  atom.  A  discussion  of  the  interplay  between  the  reduced 
mass  and  geometric  effects  for  atom- rigid  rotor  collisions  was  included  in  a 
previous  publication.^28^ 

For  kinematic  reasons,  collisions  with  the  hydrogen  atom  in  the  diatomic 
hydride  are  much  more  effective  in  transferring  precollision  center-of-mass 
kinetic  energy  into  the  diatom’s  internal  modes.  It  is  well  known  that  for  sim¬ 
ple  atom-diatomic  models,  collisions  perpendicular  to  the  diatomic  bond  activate 
only  rotation,  while  collinear  collisions  activate  only  vibration.  As  can  be  seen 
from  Eq.  (4-35),  for  impacts  on  the  hydrogen  atom  ihe  model  used  in  this  work 
predicts  that  a  specific  fraction  of  the  impulse  imparted  during  a  collision  is 
transmitted  to  the  internal  modes.  This  impulse  is  transmitted  into  the  vibra¬ 
tional  mode  for  collinear  collisions  and  into  the  rotational  mode  for  perpendic¬ 
ular  collisions,  and  is  divided  between  the  two  modes  at  inter  nediate  angles. 
Table  4-m  shows  this  fraction  of  center-of-mass  energy  f  from  Eq.  (4-36)] , 
transmitted  into  internal  modes  for  normal  collisions  on  the  diatomic  hydrogen 
atom  as  a  function  of  the  atom-diatomic  hydride  collision  pair.  These  entries, 
then,  represent  the  maximum  fraction  of  the  centor-of-mass  energy  which  can 
go  into  either  internal  mode  for  each  collision  pair...  These  results  can  be  con¬ 
trasted  with  the  fact  that  for  each  of  the  atom-rigid  rotor  pairs,  impact  points 
exist  that  allow  transmittal  of  all  the  center-of-mass  kinetic  energy  into  the 
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Table  4 -III.  Maximum  Normalized  Energy  Transfer 
to  Internal  Modes  (E^/E^j^or  Ey/EQjyj) 
in  Atom-Harmonic  Oscillator  Collisions. 


Diatomic  Atom 

He 

O 

Ar 

OH 

0.744 

0.405 

0.300 

HF 

0.730 

0.379 

0.271 

HC1 

0.691 

0.311 

0.  195 

diatomic  rotational  mode.  This  difference  between  the  dynamics  of  the  atom- 

rigid  rotor  and  the  atom-harmonic  oscillator  systems  is  mirrored  in  the  ex- 

»» 

tended  tails  of  the  PfEp/E^^)  distributions  shewn  in  Figs.  4-2  through  4-9. 

Some  uncertainty  exists  in  determining  the  proper  way  to  interpret  the 
results  of  classical  calculations  involving  the  collisional  activation  or  deactiva¬ 
tion  of  a  quantized  internal  mode  when  the  average  classical  collision  exchanges 
less  than  a  full  quantum  of  energy  with  the  mode  in  question.  Thompson  has 
discussed  this  point  recently  in  regard  to  his  Monte  Carlo  trajectory  calcula- 
tions  of  vibrational  deactivation. '  '  ”  Parker  has  shown  that  a  classical 

calculation  of  the  collinear  collisional  excitation  of  a  harmonic  oscillator  will 
yield  the  same  average  vibrational  energy  transfer  as  the  quantum  mechanical 
solution  first  proposed  by  Jackson  and  Mott  '  if  the  initial  collisional  kinetic 
energy  is  much  greater  than  the  energy  of  a  vibrational  quantum  transition  in 
the  quantum  treatment.  For  the  diatomic  hydrides  considered  in  this  work,  the 
transition  from  the  ground  to  the  first  vibrational  state  requires  on  the  order  of 
0.4  to  0.5  eV.  Thus,  unless  the  average  impulse  transmitted  along  the  diatomic 
bond  axis  becomes  significantly  larger  than  0.4  eV,  the  vibrational  energy  dis¬ 
tribution  functions,  P(E^./ECI^),  shown  in  Figs.  4-2  through  4-10,  will  be 
largely  invalid.  Since  the  great  majority  of  any  sample  of  collisions  will  be 
both  grazing  and  off  the  collinear  axis,  the  average  initial  center-of-mass 
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collisional  energies  must  exceed  the  1-  to  10-eV  range  for  which  our  calculation 
was  originally  intended  before  the  vibrational  distributions  in  Figs  4-2  through 
4-10  apply.  Table  4— II  shows  that  the  average  classical  vibrational  excitation 
varies  between  3.6%  (He  +  OH)  and  1.4%  (Ar  +  HC1)  of  the  initial  center-of- 
mass  collision  energies.  Thus,  even  in  the  most  favorable  case  (He  +  OH), 
the  initial  center -of -mass  energy  must  exceed  12  eV  before  one  vibrational 
quantity  will  be  excited  on  the  average.  Since,  on  the  average,  more  of  the 
collision  energy  goes  into  rotation  and  even  the  higher  rotational  levels  are 
much  more  closely  spaced  than  vibrational  levels,  a  similar  problem  does  not 

/oo\ 

arise  in  calculating  the  rigid-rotor  rotational  excitation. ' 

The  purpose  of  this  calculation  was  to  assess  the  appropriateness  of  the 
rigid-rotor  approximation  in  the  calculation  of  rotational  excitation  probabilities 
in  collisions  with  initial  center-of-mass  energies  of  a  few  eV.  The  results  of 
the  atom-harmonic  oscillator  calculations  show  that,  classically,  vibrational 
excitation  does  reduce  the  number  of  collisions  that  exhibit  high  degrees  of 
rotational  excitation.  However,  in  the  collision  range  between  1  and  10  eV, 
most  of  the  classically  calculated  vibrational  excitation  is  "forbidden,  "  in  the 
sense  that  it  corresponds  to  less  than  a  vibrational  quantum.  It  remains  for 
future  work  to  determine,  in  detail,  how  seriously  this  nearly  "closed"  channel 
affects  the  rotational  energy  distribution  for  this  energy  range.  In  the  mean¬ 
time,  it  seems  reasonable  to  expect  the  rigid-rotor  rotational  excitation  prob- 
abilities,  PfE^/E^^),  rather  than  the  harmonic  oscillator  collision  probabil¬ 
ities,  PfEp/E^,^),  to  more  closely  mirror  the  actual  situation  in  the  energy 
range  of  1  to  10  eV.  At  higher  center-of-mass  collision  energies,  the  harmonic 
oscillator  excitation  probabilities,  PfE^/E^^)  and  P(Ey/E^^j),  should  be 
closer  to  reality. 

The  calculations  presented  in  this  work  represent  only  a  "physical"  or 
dynamic  picture  of  inelastic  excitation  in  atom-diatomic  hydride  collisions. 

They  can  be  expected  to  be  reasonably  approximate  for  collisions  between  such 
closed-shell  species  as  Ar  and  He  with  HCl  ard  HF. 
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Although  "chemical"  effects  have  generally  been  disregarded  in  this  study, 
two  types  of  chemical  effects  can  be  expected  to  modify  the  picture  presented 
above.  The  first  is  actual  chemical  reaction,  which  should  be  most  important 
in  the  interaction  of  oxygen  atoms  with  the  hydroxyl  radical.  The  reaction 

0+  OH  —  02  +  H  (4-51) 

is  over  0.7  eV  exothermic,  while  exchange  of  the  hydrogen  atom  between  oxygen 

atoms  is,  of  course,  thermoneutral.  Reaction  (4-51)  has  a  bimolecular  rate 

“11  3  (34) 

constant  of  2. 16  x  10  cm  /sec,  apparently  independent  of  temperature. '  ' 

If  this  rate  holds  to  the  highly  nonthermal  velocities  modeled  in  this  calculation 

5 

it  would  represent  a  reaction  cross  section  (at  a  relative  velocity  of  6  x  10  cm/sec) 
2 

of  3. 6  A  .  The  approximate  hard-sphere  cross  section  for  the  models  of  O  and 

2 

OH  used  above  is  10.6  A  ,  so  that  an  appreciable  fraction  of  the  collisions  be¬ 
tween  O  and  OH  may  be  expected  to  follow  a  chemically  reactive  channel,  rather 
than  the  inelastic  activation  channel  calculated  above.  A  chemical  reaction  could 
also  be  important  in  the  nearly  thermoneutral  hydrogen  atom  exchange  between 
O  and  HC1. 

The  second  type  of  "chemical"  effect  has  been  demon°trated  in  a  number  of 
shock  tube  studies  where  nonrare-gas  atom  collisions  with  diatomics  have  shown 
an  anomolously  large  cross  section  for  vibrational  deactivation.  This 

effect  has  been  attributed  to  the  increased  "stickiness"  of  the  collisions,  as 
represented  by  the  formation  of  some  type  oi  moderately  long-lived  collision 
complex  that  can  effectively  transfer  energy  out  of  the  diatomic's  vibrational 
mode.  Nikitin  and  Umunski^*^  have  recently  speculated  that  this  stickiness  is 
due  to  "vibronic"  vibrational-to-translational  transfer  involving  excited  elec¬ 
tronic  states  of  the  atom-diatomic  complex.  Because  of  the  principle  of  micro¬ 
scopic  reversibility,  this  effect  could  also  be  important  in  the  translational-to- 
vibrational  energy  exchange  modeled  above  for  collisions  involving  O  and  OH. 


collisional  energies  must  exceed  the  1-  to  10-eV  range  for  which  our  calculation 
was  originally  intended  before  the  vibrational  distributions  in  Figs  4-2  through 
4-10  apply.  Table  4— II  shows  that  the  average  classical  vibrational  excitation 
varies  between  3.6%  (He  +  OH)  and  1.4%  (Ar  +  HC1)  of  the  initial  center-of- 
mass  collision  energies.  Thus,  even  in  the  most  favorable  case  (He  +  OH), 
the  initial  eenter-of-mass  energy  must  exceed  12  eV  before  one  vibrational 
quantity  will  be  excited  on  the  average.  Since,  on  the  average,  more  of  the 
collision  energy  goes  into  rotation  and  even  the  higher  rotational  levels  are 
much  more  closely  spaced  than  vibrational  levels,  a  similar  problem  does  not 
arise  in  calculating  the  rigid-rotor  rotational  excitation. 


The  purpose  of  this  calculation  was  to  assess  the  appropriateness  of  the 
rigid-rotor  approximation  in  the  calculation  of  rotational  excitation  probabilities 
in  collisions  with  initial  center-of-mass  energies  of  a  few  eV.  The  results  of 
the  atom-harmonic  oscillator  calculations  show  that,  classically,  vibrational 
excitation  does  reduce  the  number  of  collisions  that  exhibit  high  degrees  of 
rotational  excitation.  However,  in  the  collision  range  between  1  and  10  eV, 
most  of  the  classically  calculated  vibrational  excitation  is  "forbidden,  "  in  the 
sense  that  it  corresponds  to  less  than  a  vibrational  quantum.  It  remains  for 
future  work  to  determine,  in  detail,  how  seriously  this  nearly  "closed"  channel 
affects  the  rotational  energy  distribution  for  this  energy  range.  In  the  mean¬ 
time,  it  seems  reasonable  to  expect  the  rigid-rotor  rotational  excitation  prob- 
abilities,  P(Ej./E_,m),  rather  than  the  harmonic  oscillator  collision  probabil- 

v 

ities,  P(EJ./E.-,,,),  to  more  closely  mirror  the  actual  situation  in  the  energy 
range  of  1  to  10  eV.  At  higher  center-of-mass  collision  energies,  the  harmonic 
oscillator  excitation  probabilities,  P(E^/E^,^)  and  Pv'Ey/E^,^),  should  be 
closer  to  reality. 


The  calculations  presented  in  this  work  represent  only  a  "physical"  or 
dynamic  picture  of  inelastic  excitation  in  atom-diatomic  hydride  collisions. 
They  can  be  expected  to  be  reasonably  approximate  for  collisions  between  such 
closed-shell  species  as  Ar  and  He  with  HCl  nr«i  HF. 
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Although  "chemical"  effects  have  generally  been  disregarded  in  this  study, 
two  types  of  chemical  effects  can  be  expected  to  modify  the  picture  presented 
above.  The  first  is  actual  chemical  reaction,  which  should  be  most  important 
in  the  interaction  of  oxygen  atoms  with  the  hydroxyl  radical.  The  reaction 

O  +  OH  —  02  +  H  '4-51) 

is  over  0.7  eV  exothermic,  while  exchange  of  the  hydrogen  atom  between  oyygen 

atoms  is,  of  course,  thermoneutral.  Reaction  (4-51)  has  a  bimolecular  rate 

constant  of  2. 16  x  10  cm*3 /sec,  apparently  independent  of  temperature.  ' 

If  this  rate  holds  to  the  highly  nonthermal  velocities  modeled  in  this  calculation, 

5 

it  would  represent  a  reaction  cross  section  (at  a  relative  velocity  of  6  x  10  cm/sec) 
2 

of  3.6  A  .  The  approximate  hard-sphere  cross  section  for  the  models  of  O  and 

2 

OH  used  above  is  10.6  A  ,  so  that  an  appreciable  fraction  of  the  collisions  be¬ 
tween  O  and  OH  may  be  expected  to  follow  a  chemically  reactive  channel,  rather 
than  the  inelastic  activation  channel  calculated  above.  A  chemical  reaction  could 
also  be  important  in  the  nearly  thermoneutral  hydrogen  atom  exchange  between 
O  and  HC1. 

The  second  type  of  "chemical"  effect  has  been  demonstrated  in  a  number  of 
shock  tube  studies  where  nonrare-gas  atom  collisions  with  diatomics  have  shown 
an  anomolously  large  cross  section  for  vibrational  deactivation.  (35-39)  This 
effect  has  been  attributed  to  the  increased  "stickiness"  of  the  collisions,  as 
represented  by  the  formation  of  some  type  of  moderately  long-lived  collision 
complex  that  can  effectively  transfer  energy  out  of  the  diatomic's  vibrational 
mode.  Nikitin  and  Un.unski^^  have  recently  speculated  that  this  stickiness  is 
due  to  "vibronic"  vibrational-to-translational  transfer  involving  excited  elec¬ 
tronic  states  of  the  atom-diatomic  complex.  Because  of  the  principle  of  micro¬ 
scopic  reversibility,  this  effect  could  also  be  important  in  the  translational-to- 
vibrationai  energy  exchange  modeled  above  for  collisions  involving  O  and  OH. 
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As  noted  in  Section  3,  the  problem  of  multiple  collisions  has  also  been 
ignored  in  the  present  work.  Multiple  impacts  are  important  for  collinear  or 
nearly  collinear  collisions,  and  can  influence  the  vibrational  excitation  distri  ¬ 
bution.  Since  the  focus  of  this  work  is  on  the  rotational  excitation  distribution, 
which  is  unaffected  by  collinear  collisions  in  the  impulse  approximation,  it  was 
felt  that  ignoring  multiple  impacts  was  a  justifiable  approximation.  The 
classical  calculation  of  collinear  collisions  between  hard  sphere  atoms  and 
diatomics,  including  the  effects  of  multiple  collisions,  has  been  explored  in 
detail  by  Benson  et  al. 
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5.  PURE  ROTATIONAL  EMISSION  FROM  DIATOMIC  MOLECULES 
IN  THE  QUANTUM  AND  CLASSICAL  LIMITS 

5. 1  Introduction 

One  of  the  most  useful  concepts  in  the  construction  of  the  quantum  theory 

(42) 

is  the  correspondence  principle. '  '  It  can  be  stated  in  several  ways,  including 
the  following  way:  "Quantum  theory  must  approach  classical  theory  asymptoti¬ 
cally  in  the  limit  of  large  quantum  numbers.  Thus,  when  dealing  with  mani¬ 
festly  macroscopic  phenomena,  classical  physics  is  an  appropriate  description 
of  nature.  There  is,  therefore,  a  transition  region  where  classical  physics  is  a 
more-or-less  good  approximation  to  the  more  exact  quantum  description  of  what¬ 
ever  system  is  of  interest,  or  put  somewhat  differently,  classical  physics  is  one 
approximation  in  the  arsenal  of  all  possible  approximations  to  the  complete, 
exact  quantum  theory.  As  such,  the  quantitative  limits  on  the  region  of  validity 
of  the  classical  approximation  for  a  system  are  clearly  of  interest,  just  as  are 
the  limits  of  validity  of  any  approximation. 

Furthermore,  there  are  two  reasons  why  the  limits  of  validity  and  use¬ 
fulness  for  the  classical  theory  are  especially  interesting.  First,  a  calculation 
that  displays  the  explicit  evolution  from  quantum  to  classical  behavior  as  quantum 
numbers  increase  is  a  quantitative  illustration  of  the  correspondence  principle. 
Second,  classical  calculations  tend  to  be  more  tractable  than  the  analogous  quan¬ 
tum  calculations,  so  that  when  the  classical  theory  may  be  usefully  applied,  it  is 
usually  advantageous  so  to  do. 

The  problem  that  is  addressed  here  is  that  of  pure  rotational  radiation 
emission  from  diatomic  molecules  with  a  permanent  electric  dipole  moment. 

In  particular,  the  emission  into  a  broad  frequency  band,  which  includes  the  con¬ 
tributions  of  several  individual  transitions,  is  studied.  This  calculation  is  done 
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both  purely  quantum  mechanically  and  purely  classically;  then  the  results  of  the 
two  theories  are  quantitatively  compared  to  demonstrate  the  onset  of  classical 
behavior  as  the  relevant  quantum  numbers  (and  hence  energies)  increase. 

To  compare  the  classical  and  quantum  theories  of  pure  rotational 
emission,  it  is  necessary  to  consider  radiation  in  a  broad  wavelength  band. 

If  spectral  information,  line-by-line,  were  desired,  a  ..’'issical  approximation 
would  be  of  little  use.  Thus,  the  wavelength  band  must  enclose  a  number  of 
rotational  lines  for  straightforward  comparison  between  the  two  theories.  The 
results  of  the  calculation  show  that  for  large  quantum  numbers,  a  band  con¬ 
taining  as  few  as  approximately  five  lines  is  sufficient  for  reasonably  close  agree¬ 
ment  between  the  two  theories. 

Quantitatively ,  the  results  of  this  calculation  apply  only  to  simple  rigid 
rotors  whose  energy  levels  are  (2  f  +  l)-fold  degenerate,  where  l  is  the  princi¬ 
pal  rotational  quantum  number.  Qualitatively,  however,  these  results  shed  light, 
as  well,  on  the  classical  limit  for  the  asymmetric  rotor  with  a  permanent  electric 
dipole  moment.  Since  the  levels  of  an  asymmetric  rotor  are  nondegenerate,  and 
the  energy  levels  are  more  closely  spaced  on  the  average  than  for  the  case  of  a 
simple  rotor,  it  is  clear  that  whenever  the  classical  model  of  rotational  emission 

is  valid  for  the  diatomic  case,  the  classical  model  for  radiative  emission  from  an 

(44) 

asymmetric  rotor  will  certainly  be  valid  for  similar  quantum  numbers. 


5.2  Quantum  Mechanical  Formulation 

The  only  nonzero  Einstein  coefficients  for  purely  radiative  decay  for  a 

rigid  rotor  with  permanent  electric  dipole  moment,  p  ,  and  moment  of  inertia 

I,  are  those  in  which  the  principal  quantum  number  decreases  by  1,  and  the 

(45) 

second  quantum  number,  m,  changes  by  0  or  +  1.  ' 


A  .  %  A  tmJJ.L-m) 

fm*f-l,m  (2  f  +1)  (2  (-1) 


(5-1) 
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(5-2) 


A  =  k-  (  <«+m>  (f+™  ~  I) 

(m-(  -  1,  m-1  2  (21+1)  (21-1) 


A  ~  K  * 

*  m— 2 


(2  f+1)  (2f-l) 


where 


(41i2/x2/3cJI3) 


The  Einstein  coefficient  for  the  transition  from  any  m  value  in  state  ( 
to  all  m'  values  in  state  (  -1  is  simply  the  sum  of  the  above  three  expressions, 
and  is  independent  of  m.  Thus, 

Af  =«f4/(2f+l)  .  (5-5) 

This  means  that  the  problem  of  radiative  cascade  from  one  rotational  energy 
level  to  the  next  is  independent  of  the  quantum  number,  m. 

In  order  to  determine  the  power  radiated  by  an  ensemble  of  excited 
rigid  rotors,  the  population  of  each  excited  rotor  energy  level  must  be  found. 

If  it  is  assumed  that  at  t  =  0  there  is  an  ensemble  of  rotors  all  excited  into  the 
rotational  state  with  principal  rotational  quantum  number  f  ,  then  the  equation 
governing  the  average  population  of  the  f  level  is  given  by  the  following 
expression,  in  the  absence  of  collisions: 


-  =  “Af  Nf  (t) 
o  o 


(5-6) 


where  N.  (t)  is  the  average  population  of  the  l  level,  regardless  of  the  quantum 
lo  0 

number  m.  If  the  normalized  boundary  condition  of  NfQ(t=0)  =  1  is  imposed,  the 
solution  of  this  simple  equation  is 


-A,  t 


Nf0  (t)  =  e 


(5-7) 


In  general,  the  population  of  the  f  -n  level  obeys  the  following  equation: 


dN 


f  -n 
o 


dt 


‘A1  -n  N/  -n 
o  o 


A,  N, 

o-n+1  o - n  1 


(5-8) 


Because  of  the  selection  rule,  f'  =  f  -  1  for  radiative  decay,  only  the  population 
of  the  level  one  step  higher  than  the  level  of  interest  can  contribute  to  the  popu¬ 
lation.  Thus,  the  only  two  nonzero  terms  on  the  right  side  of  the  above  equation 
have  been  included. 

The  above  set  of  equations,  along  with  the  equation  and  the  boundary 
condition  at  t  =  0,  is  easily  solved.  For  1 1  fQ,  the  general  solution  is 


Nf(t) 


2 


#*=f 


-Af,  t 
e 


Af  Af  -1  •**  Af+i 
v  o  o 


fo 

n  (a 

f"  =  t 
p”  /n 


p " 


(5-9) 


The  above  equation  gives  the  average  population  of  the  various  rotational  energy 

levels  for  an  ensemble  of  rotors  with  the  normalized  boundary  condition  that 

N  (t=0)  =  1  and  Nf  ,  f  (t=0)  0. 

fo  ^  o 
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The  energy  given  up  by  a  rigid  rotor  when  it  undergoes  a  radiative  transi¬ 
tion  from  state  (  to  state  «  -1  is 

Ef  -  Ef  l  -  (h2/2I)rf(f+l)  -  (f-1)  ((-1  +  1)]  =  (h2/I)  /  (5-10) 

th 

The  contribution  to  the  average  radiated  power  by  the  f  level  Pf  is,  therefore, 
given  by 

Pf  =  N{(t)  (E,-Ef  l)  Af  .  (5-11) 

To  compute  the  total  average  emitted  power,  it  is  now  only  necessary  to  sum 
up  the  contribution  from  each  individual  transition.  The  average  total  power 
emitted  by  an  ensemble  of  particles,  all  of  which  were  in  the  state  f  at  t  =  0,  is 


The  above  expression  is  the  total  radiated  power.  If  only  the  power  emitted 
into  some  frequency  band  is  of  interest,  the  sum  overf  must  be  modified  to  include 
only  those  transitions  that  contribute  to  the  band.  If  the  lower  and  upper  limits 
for  the  angular  frequency  for  the  band  are  designated  and  respectively, 
then, 

S  (E(  -  Em)/H<u>2  .  (5-13) 


1-97 


This  is  an  implicit  relation  for  the  maximum  and  minimum  allowed  f 's  in 
the  first  sum  over  ( in  Eq.  (5-12).  The  power  in  a  given  band  will  simply  be 
denoted  by  P  (  rQ  ,  t;  ojj  , 

Equation  (5-12)  gives  the  power  radiated  as  a  function  of  time.  Therefore, 
to  find  the  total  energy  radiated  into  the  band,  P  must  be  integrated  with  respect 
to  time  from  t  =  0  to  t  =  ®  .  Thus, 


EN  ( cj2) 


dt  P  (  fQ,  t;  <ur  tu2) 


(5-14) 


and 


EN  <fo;  wl*  U)2>  = 


f  -1 


eoit/i-Uj)  e(w2-h  f  /i) 


+  0(hf  /I-W  )  6(01-1} e  /I)  f 
'  o  1  2  o  o 


(5-15) 


where  $  (x) 


0  for  x  <  0 
1  for  x  >  0 


(5-16) 


The  two-step  functions  in  Eq.  (5-15)  arise  because  of  the  band  limits  on  the  first 
f  sum  (i.e.,  as  a  result  of  Eq.  (5-13). 

Rewriting  the  sum  of  the  products  of  the  pair-wise  differences  of  the  A{  's 
in  Eq.  (5-15)  in  terms  of  a  Vandermonde  determinant,  ^^then  putting  all  the 
terms  in  the  f  sum  over  this  common  denominator,  the  total  band  energy',  EN, 
simplifies  considerably. 
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( (o\(jJV0J2 )  --=  (h2/l>  f  Bfat/I-Uij  9(UJ?  -t»f/l)  .  (5-17) 


This  expression  simply  means  that  each  level  in  the  band  pass  contributes  an 
energy  of  (h  f/l)  to  the  total  energy,  which  could  have  been  written  in  the 
beginning-  The  integration  of  the  power  to  get  this  simple  expression  for  the 
energy  is,  therefore,  a  check  on  the  validity  of  the  power  formula,  Eq.  (5-12). 

Equations  (5-12)  and  (5-17)  are  the  quantum  mechanical  expressions 
for  the  power  and  energy  emitted  by  a  rigid  rotor  with  a  permanent  electric  dipole 
moment.  The  power  emitted  by  such  a  rotor  may  also  be  calculated  by  using 
strictly  classical  physics.  It  is  this  classical  version  of  the  calculation  to  which 
attention  is  now'  directed. 


5.3  Classical  Formulation 

A  classical  rigid  rotor  with  a  permanent  electric  dipole  moment  radiates 

as  it  rotates,  since  the  two  charges  making  up  the  dipole  are  being  accelerated 

because  of  their  circular  motion.  The  dipole  is  radiating  energy  away,  so  that 

(44  47) 

there  is  an  effective  torque  acting  on  it,  the  radiative  damping  torque,  ry  ’ 


r  =  (2/3c  )  x  p 


(5-18) 


where  /x  is  the  electric  dipole  moment,  and  the  dots  signify  a  time  derivative. 
For  pure  rotation 


=  fcl  XU  X  u  X  U  , 

*»•  ** 

where  u;  is  the  rotor's  angular  velocity. 


(5-19) 
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This  classical  model  for  radiative  damping  and  emission  now  reduces  to 
solving  the  equation  of  motion  for  <u  ,  using  the  torque,  r.  For  a  simple 
rigid  rotor, 

I  cu  =  r  ,  (5-20) 

M 

or 

w  =  (2^2/3c3D  a,3  ,  (5-21) 

The  solution  of  this  equation  is  straightforward,  and  leads  to  the  following 
expressions  when  using  the  initial  condition  that  (u  (t=0)  =  <u 

2  1/2 

cuo/(l  +  6  uj^t)  ,  (5-22) 

where 

2  .  3 

6  -  4 ix  /3c  I  .  (5-23) 

In  order  to  find  the  power  radiated  by  this  rotating  dipole,  it  is  necessary 
only  to  use  these  expressions  for  (u  in  Larmor's  power  formula: 

P  =  (2/3c3)  |  |  2  ,  (5-24) 

or,  upon  using  Eq.  (5-22), 

P  =  (16/2)  (Uq4/  (1  +6(uq2  t)2  .  (5-25) 

This  is  the  classical  analog  to  Eq.  (5-12)  for  the  power  emitted  from  a  quantum 
mechanical  rigid  rotor  with  a  permanent  electric  dipole  moment. 

The  above  equation  is  the  total  average  power  emitted  by  the  rotating  dipole. 
If  one  is  interested  in  the  power  emitted  into  some  frequency  band  from  <u  to 
(Uj  with,  for  example,  uig  >  oj  j  ,  then  this  expression  must  be  modified  to 
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include  contributions  between  only  the  proper  frequencies.  Thus, 


P(wru»2)  =  [(16/2)  wq4/  (1  +  26  o»o2t)] 

x  [e  («2  -<*>)  -  Wj)]  .  (5-26) 

The  total  energy  radiated  during  a  complete  band  pass,  EN,  by  a  classical 
rotor  is  simply  the  difference  in  energy  of  the  rotor  between  the  state  when  it  is 
at  the  top  of  the  band  and  the  state  when  it  is  at  the  bottom  of  the  band. 


EN  =  I  (1/2)  I(  u>22  -  u2)  for  >  u>2 


(1/2)  I(  u)q2  -  Wj2)  for  u2  <  u>Q  <  w. 


(5-27) 


for  u  <  u, 
o  1 


In  the  completely  classical  limit,  the  power  described  by  Eqs.  (5-12)  and 
(5-26)  should  become  identical.  In  fact,  if  a  classical  approximation  is  made  in 
Eq.  (5-8),  the  classical  result  in  Eq.  (5-22)  may  be  obtained  directly.  For  a 
state,  f  ,  Eq.  (5-8)  may  be  rewritten  as 


2  f  +  1 


(2  (f-  l)  +1)  1  (-1 


(5-28) 


having  used  the  expression  for  the  Einstein  coefficient,  Eq.  (5-5).  In  order  to 
approach  the  classical  limit  in  the  equation,  l  must  become  a  continuous  variable 
so  that  N  becomes  a  function  of  the  two  variables  t  and  £ .  Also,  the  difference 
between  the  f  ^  and  the  (f  +  1)^  terms  on  the  right  side  of  Eq.  (5-28)  must 
become  a  partial  derivative.  Finally,  by  recognizing  that 


2  f  +  1  =  2f 


(5-29) 
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in  the  classical  limit,  Eq.  (5-28)  becomes 


1  _JL 

2  3f 


[f3  N  (f  ,  t)] 


(5-30) 


f  <3 


f*  N 


(5-31) 


The  characteristic  equations  for  this  partial  differential  equation  are 


(5-32) 


d£  _ 
ds 


-I  '3 


(5-33) 


The  solution  of  these  equations  is  straightforward. 


(5-34) 


t  -  s/k 


(5-35) 


f2  -  (S.^V1 


(5-36) 


Having  used  the  boundary  condition  that  when  s  is  zero  (and  hence,  when  t  is  zero), 
f  is  f  .  From  Eq,  (5-10),  the  relation  between  the  frequency  cj  and  the  quantum 
number  f  in  the  classical  limit  is  clearly  given  by 


w  *  tit  /l. 


(5-37) 
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Thus,  using  Eq.  (5-35)  along  with  the  definitions  of  k  and  6  ,  Eq.  (5-36)  becomes 
U)2  =  WQ2/  (1  +  6  u>o2t)  ,  (5-38) 

which  is  identical  to  the  purely  classical  result,  Eq.  (5-22). 

5.4  Results 

In  order  to  examine  the  onset  of  the  classical  limit  for  the  rotational  radi¬ 
ation  from  a  diatomic  molecule,  the  power  as  calculated  from  both  quantum 
(Eq.  (5-12))  and  class ical^Eq.  (5-26))  theories  has  been  computed  for  a  variety 
of  initial  quantum  levels  and  band  passes.  The  curves  for  radiated  power  as  a 
function  of  time  are  shown  in  Figs.  5-1  through  5-4.  In  each  case  the  prediction 
of  the  quantum  theory  is  represented  by  the  solid  lire,  and  the  prediction  of  the 
classical  theory  is  represented  by  the  dashed  line.  The  classical  curve  has 
been  drawn  so  that  the  area  under  the  curve  (i.e. ,  the  total  energy  emitted  in  a 
band  pass)  is  the  same  as  that  for  the  quantum  curve.  This  means  that  the  energy 
of  the  rotating  molecule  as  it  decays  into  the  band,  classically,  is  the  same  as 
the  highest  quantum  energy  state  included  in  the  band;  and  for  the  decay  out  of  the 
band,  the  classical  energy  is  equal  to  the  lowest  quantum  state  included  in  the  band. 

The  expressions  for  the  radiated  power  are  applicable  to  any  diatomic 
molecule  with  a  permanent  electric  dipole  moment,  but  in  order  to  be  concrete, 
all  the  results  have  been  plotted  for  the  hydrogen  flouride  molecule.  The  results 
will  be  qualitatively  the  same  for  other  molecules,  but  since  the  moment  of  inertia 
and  electric  dipole  moment  will  vary  from  case  to  case,  the  time  scale  for  decay 
will  be  different  for  each  case.  It  is  a  reasonably  straightforward  matter  to  apply 
the  results  to  any  other  molecule  of  interest. 

There  are  a  number  of  interesting  features  shown  in  Figs.  5-1  through  5-4. 
First  of  all,  as  the  relevant  quantum  numbers  increase,  the  agreement  between 
the  quantum  and  classical  theories  improves,  as  one  would  expect  on  the  basis  of 
the  correspondence  principle.  In  fact,  the  results  plotted  in  Fig.  5-3  for  an 
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Figure  5.1-  Radiated  power  as  a  function  of  time  for  the  molecule  HF  with  initial  rotational  quantum 
number  %  =10,  and  for  a  band  pass  including  contributions  from  the  6th  through  the  4th 


Figure  5.2  -  Radiated  power  as  a  function  of  time  for  the  molecule  HF  with  initial  rotational  quantum 
number  9  =15,  and  for  a  band  pass  including  contributions  from  the  10th  through  the  5th 


Figure  5.4  -  Radiated  power  as  a  function  of  time  for  the  molecule  HF  with  initial  rotational  quantum 
number  9  =20,  and  for  a  band  pass  including  contributions  from  the  20th  through  the  lot 

’  r\  * 


initial  rotational  quantum  number  of  30,  and  a  band  enclosing  contributions  from 
the  25th  through  the  10th  levels,  sho’"  close  agreement  between  the  two  theories 
for  almost  all  times  of  interest. 

In  Figs.  5-1  through  5-3,  the  initial  rotational  state  emits  at  a  frequency 
too  high  to  contribute  to  the  band  power,  so  that  there  is  some  nonzero  time  inter¬ 
val  when  there  is  little  or  no  radiation.  Classically,  the  onset  of  radiation  is 
abrupt,  but  quantum  mechanically,  it  is  gradual,  and,  in  fact,  is  nonzero  before 
the  beginning  of  the  classical  radiating  region.  This  precursor  is  a  manifesta¬ 
tion  of  the  fact  that  there  is  always  a  nonzero,  albeit  sometimes  quite  small, 
probability  for  a  very  rapid  quantum  transition. 

At  the  onset  of  the  classical  radiation,  the  classically  computed  power  is 
larger  than  the  quantum  mechanically  computed  power  by  roughly  a  factor  of  1.5 
or  2.  Tlr  ,  "overshoot"  remains  at  higher  quantum  levels,  but  tends  to  become 
increasingly  narrow,  so  that  this  initial  classical  peak  becomes  a  spike  of  infini¬ 
tesimal  area. 

Finally,  just  as  the  classical  curve  begins  abruptly,  it  ends  abruptly  as 
the  frequency  decays  to  a  value  outside  the  band  limits.  The  quantum  curve 
decays  continuously  across  this  region,  and  the  decay  is  asymptotically  described 
by  an  exponential  function  with  a  lifetime  equal  to  the  lifetime  of  the  lowest 
energy  state  that  contributes  to  the  band. 

As  Fig.  5-1  through  5-3  show,  the  main  difference  between  the  classi¬ 
cal  power  and  the  quantum  power  for  quantum  numbers  above  approximately 
10  occurs  in  the  precursor  and  overshoot  at  the  beginning  of  the  classical  radi¬ 
ating  region.  Figure  5-4  shows  the  case  where  the  initial  state  is  able  to  con¬ 
tribute  to  the  band  power,  so  that  there  is  neither  a  precursor  nor  an  overshoot. 
The  agreement  between  the  two  curves  is  within  a  few  percent  everywhere  in  the 
classical  radiating  region. 
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A  number  of  conclusions  can  be  drawn  from  the  foregoing  discussion. 

For  rotational  quantum  numbers  greater  than  15  to  20,  and  for  bands  sufficiently 
wide  to  include  the  contributions  of  6  or  more  levels,  the  classical  theory  very 
closely  reproduces  the  quantum  mechanical  prediction  of  the  radiated  power  as  a 
function  of  time.  This  is  especially  true  in  those  cases  in  which  the  initial  state 
can  contribute  to  the  power  radiated  into  the  band.  Any  distribution  of  initial 
radiating  states  (e.g. ,  a  thermal  distribution)  would  remove  the  abrupt  beginning 
and  end  of  the  classical  radiating  region.  Hence,  the  sharp  corners  of  the  classi¬ 
cal  curve  would  be  smoothed,  bringing  the  predictions  of  the  two  theories  into 
considerably  closer  agreement. 

This  work  is  for  the  case  of  a  simple  or  diatomic  rotor  whose  energy 
levels  are  (2  l  +l)-fold  degenerate.  For  an  asymmetric  rotor  whose  levels  are 
not  degenerate,  the  classical  theory  would  probably  be  adequate  at  even  lower 
quantum  numbers  or  energies,  and  for  bands  somewhat  narrower  than  those  con¬ 
sidered  for  the  diatomic  rotor.  The  classical  theory  for  an  asymmetric  rotor 
rloes  not  predict  sharp  cutoffs  for  the  power thus,  there  would  be  relatively 
few  precursor  or  overshoot  problems,  as  described  above. 
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